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2.1  Introduction 

The  availability  of  relatively  inexpensive  and  powerful  computing  technology 
has  profoundly  changed  the  way  in  which  modern  scientific  research  is  con¬ 
ducted  in  numerous  fields.  At  a  rather  trivial  but  highly  relevant  practical 
level,  the  rapid  increase  in  computer  power  has  considerably  sped  up  the 
pace  of  theoretical  schemes  and  approaches  for  simulating  bulk  material  at 
the  atomic  level.  The  results  of  the  simulations  have  been  invaluable  in  the 
guidance  of  experiments  and  for  providing  insight  into  system  behavior,  par¬ 
ticularly  under  extreme  conditions  of  temperature  and  pressure.  In  this  re¬ 
spect,  a  major  role  has  been  played  by  molecular  dynamics  (MD),  which  since 
the  pioneering  efforts  of  the  1960s  [1,2]  has  developed  into  a  mature  and 
active  discipline  that  has  been  used  as  a  means  of  simulating  and  understand¬ 
ing  the  properties  of  real  systems.  More  recently,  major  progress  has  been 
achieved  in  the  development  of  ab  initio  or  “first  principles”  MD,  in  which 
the  potential  energy  and  interatomic  forces  are  derived  from  accurate  quan¬ 
tum  mechanical  electronic  structure  calculations  that  are  performed  as  the 
simulation  proceeds  [3-5].  This  has  greatly  improved  the  predictive  power  of 
the  simulation  and  opens  the  way  for  the  reliable  simulation  of  processes  in 
which  chemical  bonds  are  formed  and  broken.  The  formation  and  breaking 
of  chemical  bonds  is  simulated  with  great  difficulty  in  MD  based  on  classical 
mechanics  and  empirically  derived  reactive  potentials  [6-9]. 

The  development  of  ab  initio  MD,  along  with  the  rapid  development  of 
computer  power,  has  inspired  its  application  to  studies  of  energetic  materi¬ 
als  and  detonation  processes.  There  are  several  difficult  problems  currently 
encountered  in  the  field  of  molecular  condensed  matter  detonation.  Is  the  ini¬ 
tiation  of  detonation  in  existing  molecular  explosives  controlled  by  the  ther¬ 
mal  decomposition  via  shock  temperature  or  by  initial  molecular  collisions 
occurring  within  the  shock  front?  [6, 10]  Could  the  initial  molecular  collision 
and  subsequent  bond  energy  release  be  made  the  main  mechanisms  for  the 
initiation  and  propagation  of  detonation  in  a  new  generation  of  molecular 
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explosives,  whose  detonation  velocity  could  be  several  times  that  of  current 
explosives?  Would  pressure  dissociation  or  momentum-induced  decomposi¬ 
tion  at  low  temperatures  be  another  possible  mechanism  for  detonation  to 
generate  high-speed  solid  particles  products?  Experimental  support  to  an¬ 
swer  these  questions  must  be  derived  under  extreme  conditions  of  pressure 
and  temperature  from  observations  inside  the  shock  front,  thereby  requiring 
measurements  on  the  time  scales  of  10_2-10°ps.  Molecular  simulation  has 
offered  an  alternative  avenue  towards  gaining  insights  into  these  challenging 
problems  and  guidance  to  experimental  research  with  energetic  materials  and 
detonation  processes.  High-energy-density  materials  must  feature  metastabil¬ 
ity  and  a  large  energy  content  that  mostly  originates  in  transformations  of 
the  molecular,  atomic  and  electronic  structures.  Successful  synthesis  of  these 
materials  and  control  of  the  energy  release  from  the  structural  bonds  strongly 
rely  on  the  understanding  of  the  chemical  processes  and  physics  of  structural 
transformations  at  the  atomic  level. 

In  the  following  we  review  computations  that  have  recently  been  carried 
out  in  our  group  on  high  energy  density  materials.  Quantum  mechanics- based, 
density  functional  theory  (DFT)  [11]  methods  are  used  to  provide  atomic- 
level  insight  into  the  physical  and  chemical  properties  and  processes  which 
occur  in  these  materials.  Compared  to  classical  molecular  scale  simulation 
methods  that  are  based  on  empirically  derived  potentials,  the  simulations 
presented  here  allow  for  the  simulation  of  complicated  covalent  bond  breaking 
and  formation,  such  as  in  the  case  of  molecular  decomposition. 

The  cases  discussed  include  the  reaction  of  nitromethane  under  binary 
collision  conditions,  in  the  absence  and  presence  of  spectator  molecules  which 
participate  in  secondary  collisions  in  the  system.  Additionally,  the  reaction  of 
liquid  nitromethane  under  high  pressure  conditions  is  studied.  Our  work  on  the 
solid-state  phases  of  covalent  nitrogen  is  also  summarized.  There  are  several 
areas  where  computations  can  contribute  in  these  studies.  First,  direct  infor¬ 
mation  about  optimized  crystallographic  lattice  parameters  and  the  geometri¬ 
cal  parameters  of  systems  can  be  determined.  This  can  be  done  for  a  variety  of 
uncompressed  lattices  or  when  external  isotropic  or  anisotropic  compression 
is  applied  to  the  crystal.  Besides  the  structural,  elastic,  and  phonon-modes 
parameters,  other  important  energetic  and  electronic  properties  can  be  eval¬ 
uated.  Among  the  list  of  electronic  properties  some  representative  examples 
are  the  band  structure  and  the  total  or  partial  density  of  states.  Further¬ 
more,  additional  insight  can  be  obtained  from  population  analyses  of  charge 
distribution,  bond  order,  and  electron  and  spin  density  maps. 

2.1.1  Introduction  to  First  Principles  Molecular  Simulation 

When  studying  processes  at  the  molecular  level,  the  most  fundamental  prop¬ 
erty  of  interest  is  the  molecular  potential  energy  surface  or  simply  the  poten¬ 
tial  energy  surface  (PES).  The  potential  energy  surface  of  a  system  relates 
the  total  energy  of  a  molecular  system  with  the  geometric  arrangement  of 
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Fig.  2.1.  Generic  potential  energy  surface  for  a  chemical  reaction.  The  reaction 
coordinate  represents  a  combination  of  nuclear  coordinates  which  links  the  reactant 
molecules  on  the  potential  energy  surface  to  the  products 


the  atoms  that  make  up  the  system.  Since  each  atomic  nucleus  has  three 
spatial  degrees  of  freedom  associated  (z:  y,  z  in  Cartesian  coordinate  space) 
the  potential  energy  surface  is  a  3N  dimensional  hypersurface  where  N  is  the 
number  of  atoms  that  make  up  the  system.  A  simple  one-dimensional  potential 
energy  surface  of  a  generic  chemical  reaction  is  depicted  in  Fig.  2.1. 

In  principle,  if  the  potential  energy  surface  for  a  given  reaction  can  be  cal¬ 
culated,  important  quantities  such  as  the  reaction  enthalpy,  which  is  related 
to  the  energy  content  of  a  material,  the  reaction  barriers,  the  reaction  mech¬ 
anism,  and  the  structure  of  a  material,  among  other  quantities,  can  be  deter¬ 
mined.  Thus,  a  detailed  understanding  of  a  reactive  process  can  be  achieved 
via  computer  modeling  of  the  potential  energy  surface  that  in  turn  provides 
an  avenue  for  rationally  studying  reactions  and  predicting  new  products  at 
different  temperature  and  high  pressure  conditions. 

A  wide  array  of  chemical  simulation  methods  exist  and  have  been  used  for 
decades,  as  outlined  in  Fig.  2.2.  For  processes  which  do  not  involve  the  forma¬ 
tion  or  breaking  of  chemical  bonds,  the  motion  of  reactants  on  the  potential 
energy  surface  can  be  adequately  described  by  classical  potentials.  These  are 
also  referred  to  as  empirical  potentials,  molecular  mechanics  methods,  force 
held  methods  and  ball  and  spring  potentials.  For  these  processes,  the  inter¬ 
ned  intramolecular  forces  between  the  atoms  can  be  described  loosely  as  forces 
among  classical  springs  (harmonic  or  anliarmonic).  The  force  constants  among 
these  “springs”  are  determined  empirically  by  comparing  the  predictions  of 
the  simulations  of  thermodynamic  properties  with  experimental  results,  or 
higher  level  ab  initio  calculations.  Thus  for  processes  such  as  melting,  evapora¬ 
tion,  and  non-reactive  shockwaves,  classical  potentials  may  be  used.  Classical 
potentials  have  the  advantage  that  they  are  relatively  inexpensive  to  calcu¬ 
late  and,  as  a  result,  simulations  with  thousands  to  millions  of  atoms  can  be 
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Fig.  2.2.  Spectrum  of  methods  to  calculate  the  potential  energy  surface.  Moving 
toward  the  left,  methods  tend  to  be  faster,  allowing  for  more  atoms  to  be  simulated. 
Moving  toward  the  right,  methods  are  more  transferable,  and  often  more  accurate 


performed.  The  main  disadvantages  of  classical  potentials  are  that  they  can 
not  be  generally  used  to  simulate  bond  breaking/ forming  processes,  and  they 
are  limited  in  their  transferability.  The  latter  caveat  means  that  the  empirical 
constants  in  the  potentials  are  chosen  to  reproduce  a  certain  set  of  properties 
of  a  limited  given  class  of  materials,  and  as  a  result  these  potentials  may  be 
of  limited  value  when  simulating  materials  and  properties  different  from  the 
training  set.  For  simulating  processes  which  involve  the  formation  and  break¬ 
ing  of  covalent  bonds,  the  use  of  classical  potentials  is  limited  in  value,  and  in 
these  cases,  quantum  mechanical  methods  are  most  often  used.  We  note  that 
there  are  some  ‘empirical5  potentials  that  do  allow  for  covalent  bond  breaking 
and  formation  [12-14],  but  these  are  still  very  limited  to  the  class  of  materials 
for  which  the  potentials  are  parameterized. 

With  quantum  mechanical  based  methods,  there  is  an  attempt  to  solve 
for  the  Schrodinger  equation  or  equivalent  equations,  which  provide  the  elec¬ 
tronic  structure  of  the  system  either  in  the  form  of  a  wave  function  or  electron 
density.  The  detailed  electronic  structure  allows  one  to  evaluate  the  energy 
of  the  system,  the  forces  exerted  on  the  nuclear  centers,  and  therefore  the 
all-important  potential  energy  surface.  There  exists  a.  hierarchy  of  quantum 
mechanics  methods  used  in  molecular  simulation,  ranging  from  fast,  empir¬ 
ical  methods  to  more  computationally  demanding  ‘first  principles1  methods. 
The  choice  of  method  to  be  used  in  a  simulation  generally  involves  striking 
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a  balance  between  accuracy  and  computational  cost  and  typically,  these  two 
features  of  a  method  run  contrary  to  one  another.  That  is,  there  is  usually  a 
high  computational  price  to  pay  for  high  accuracy. 

In  this  work,  we  will  restrict  ourselves  to  simulations  using  ‘first  principles’ 
quantum  mechanical  methods.  First  principles  methods  sometimes  referred  to 
as  ab  initio  methods,  use  little  or  no  empirical  data  to  construct  the  electronic 
structure.  As  a  result,  they  have  a  high  degree  of  transferability  since  they  arc 
not  constrained  by  a  particular  training  set.  This  makes  first  principles  quan¬ 
tum  mechanical  methods  well  suited  to  simulating  materials  under  extreme 
conditions  where  available  experimental  data  can  be  limited  or  non-existent. 
The  first  principles  methods  of  choice  for  simulating  materials  and  molecules 
are  density  functional  theory  (DFT)  methods.  DFT  strikes  a  favorable  balance 
between  accuracy  and  computational  efficiency  that  has  led  to  its  widespread 
use  in  the  molecular  and  material  sciences.  DFT  calculations  can  routinely 
be  performed  on  systems  containing  hundreds  or  thousands  of  atoms  with 
modern  computers. 

2.1.2  Density  Functional  Theory 

The  basic  premise  for  DFT  is  that  all  ground  state  properties  of  a  chemical 
system  are  uniquely  determined  by  the  total  electron  density  rather  than  the 
full  multielectron  wavefunction.  This  is  known  as  the  first  Hohenberg-Kohn 
theorem  [15].  This  theorem  implies  that  the  electronic  energy  of  a  molecu¬ 
lar  system  is  a  functional  of  the  electron  density.  A  functional  is  simply  a 
function  of  a  function,  for  example,  the  value  of  a  variable  /  has  a  functional 
relationship  with  integrable  functions  /(:/;),  shown  as  /[/]  for  the  following 
case,  1  [/]  =  f{x)dx .  The  value  of  I  depends  on  the  function  f{x).  In  this 

case,  the  energy,  E[p],  is  a  functional  of  the  electron  density,  p(r),  which  is 
in  turn  a  function  of  the  three  spatial  coordinates  (r  or  x}  y,  z  in  Cartesian 
space).  From  the  electronic  energy  of  a  system,  the  total  energy  of  the  sys¬ 
tem  can  be  easily  determined  for  a  given  geometric  arrangement  of  the  atoms 
that  make  up  the  system,  and,  therefore,  the  potential  energy  surface  can  be 
evaluated.  The  electronic  energy  functional  is  typically  divided  into  compo¬ 
nents  that  account  for  the  electron  kinetic  energy,  electron-nuclear  interaction 
potential  energy  and  the  electron-electron  interaction  potential  energy: 

E[  p]  =  T[p]  +  £?Ne[p]  +  1  (2-1) 

kinetic  energy  nuclear  — electron  electron  —  electron 

of  the  electrons  interaction  energy  interaction  energy 

whore  p  is  the  electron  density  and  the  square  brackets  indicate  a  functional 
dependence  in  which  the  variable  p  is  a  function  of  other  variables,  in  this 
case  the  spatial  coordinates. 
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In  1965,  Kohn  and  Sham  recognized  that  the  utility  of  (2.1)  was  severely 
limited  by  the  lack  of  a  kinetic  energy  functional,  T[p]  [16].  In  fact,  previous 
attempts  to  use  DFT  with  approximate  kinetic  energy  functionals  [17,18], 
predicted  that  chemical  bonding  should  not  occur  at  all  [19].  To  tackle  this 
situation,  Kohn  and  Sham  ingeniously  introduced  a  ‘fictitious5  non-interacting 
reference  system  (analogous  to  Hartree-Fock  orbitals)  built  from  a  set  of  single 
electron  wave  functions,  or  Kohn-Sham  orbitals,  {^(r)},  that  sum  to  give  the 
real  electron  density  p, 

P(r)  =  £hMr)|2  (2.2) 

i 

where  the  summation  index,  i,  runs  over  occupied  orbitals  (loosely  one  orbital 
for  each  electron  in  the  system).  The  kinetic  energy  of  the  non-interacting 
system  can  then  be  calculated  as: 

ro  =  E  /^(r)(-^V2)  v;t;(r)dr,  (2.3) 

Using  (2.3)  it  is  possible  to  rewrite  (2.1)  as: 

£[p]  =  To  +  £Ne[p]  +  £CG[p]-  (2.4) 

The  first  term  on  the  right  hand  side  of  (2.4)  is  the  kinetic  energy  of  the 
fictitious  electronic  system.  It  should  be  noted  that  (2.4)  is  formally  and  con¬ 
ceptually  different  than  (2.1).  In  (2.4),  the  kinetic  energy  term,  T0  refers  to 
the  kinetic  energy  of  a  non-interacting  reference  system,  whereas  in  (2.1)  the 
kinetic  energy  term  refers  to  the  true  interacting  system.  Using  the  Born- 
Oppenheimer  approximation,  through  which  all  nuclei  are  treated  as  fixed 
point  charges,  the  nuclear-electron  term  is  given  by: 

E»-lP]  =  -£/|frTj*.  (2-5) 

where  the  index  I  runs  over  all  nuclei,  Z\  is  the  nuclear  charge  of  atom  I, 
and  Rj  are  the  nuclear  coordinates  of  atom  I.  This  x*elationship  describes  the 
Coulombic  attraction  between  the  electrons  and  the  nuclei. 

The  last  term  in  (2.4),  #ee[p]>  accounts  for  the  interactions  between 
electrons.  In  Kohn-Sham  DFT  (KS-DFT),  this  term  is  separated  into  two 
components: 

em=\h~ (2.6) 

The  first  term  on  the  right  hand  side  of  (2.6)  represents  the  classical  Coulombic 
electronic  repulsion  of  an  electron  moving  in  the  average  electric  field  of  the 
other  electrons  in  the  system.  The  second  term,  /?xc,  is  the  so-called  exchange 
correlation  energy,  which  will  be  discussed  further  below. 
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When  (2.5)  and  (2.6)  are  inserted  into  (2.4),  one  arrives  at  the  following 
expressing  for  the  electronic  energy  of  the  system: 


£[p]  =  j  mo  (-^v2)  mo* 

$*$**+* «w-  m 


This  equation  forms  the  basis  of  Kohn-Sham  DFT.  In  order  to  solve  this 
equation,  one  needs  to  find  the  Kohn-Sham  wave  functions  &s  that  minimize 
the  energy  of  the  system,  subject  to  the  constraints  that  the  total  number 
of  electrons  in  the  system  is  conserved  and  the  orbitals  used  to  construct 
the  fictitious  electronic  system  are  orthogonal.  This  involves  self-consistently 
solving  a  set  of  one-electron  equations.  The  details  of  this  procedure,  although 
straightforward,  are  beyond  the  scope  of  this  discussion  and  will  not  be  de¬ 
scribed  here;  however,  the  interested  reader  is  directed  to  Chap.  7  of  [11]. 

The  first  three  terms  on  the  right  hand  side  of  (2.7)  are  well-defined  for 
a  given  set  of  orbitals.  The  accuracy  of  KS-DFT  is  thus  dependent  upon 
the  expression  for  the  last  term,  Exc ,  which  serves  as  a  depository  for  con¬ 
tributions  to  the  electronic  energy  that  are  not  adequately  described  by  the 
other  terms  in  (2.7).  These  energy  contributions  include  those  arising  from 
the  correlated  motion  of  electrons  due  to  instantaneous  electron-electron  in¬ 
teractions  (correlation),  the  correlated  motion  of  electrons  due  to  the  Pauli 
exclusion  principle  (exchange),  a  correction  for  the  difference  between  the 
kinetic  energy  of  the  fictitious  electronic  system  of  non-interacting  orbitals 
and  that  of  the  true  system,  and  a  correction  for  the  interaction  of  the  elec¬ 
tron  with  itself.  The  specific  details  of  these  energy  contributions  will  not  be 
considered  in  greater  detail  here. 

The  exchange-correlation  energy,  Exc ,  is  evaluated  with  the  use  of  a  func¬ 
tional,  which  is  called  the  exchange-correlation  (XC)  functional.  The  form  of 
this  functional  is  unknown;  however,  if  it  was  determined  exactly,  (2.7)  would 
provide  the  exact  ground  state  energy  for  the  system.  Although  the  true  XC 
functional  remains  unknown,  numerous  attempts  have  been  made  to  develop 
approximate  XC  functionals  (e.g.,  see  Chap.  6  of  reference  [11]).  In  general, 
approximate  XC  functionals  are  developed  through  either  a  consideration  of 
the  fundamental  physics  pertaining  to  electron-electron  interactions  or  fitting 
parameters  in  various  functional  forms  to  experimental  data.  Although  the 
former  approach  is  more  elegant,  the  latter  often  produces  functionals  that 
are  very  accurate  for  chemical  applications. 

The  simplest  XC  functionals  are  based  solely  upon  the  value  of  the  elec¬ 
tron  density  [20,21].  This  approach  is  termed  the  local  density  approximation 
(LDA)  and  is  generally  too  inaccurate  for  use  in  most  chemical  applications. 
Significant  improvements  in  accuracy  are  achieved  if  the  XC  functional  de¬ 
pends  on  the  value  and  gradient  of  the  electron  density  [22].  This  is  termed 
the  generalized  gradient  approximation  (GCA).  GGA  XC  functionals  are  used 
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commonly  in  many  chemical  applications.  More  sophisticated  functionals  can 
be  obtained  by  incorporating  higher-order  derivatives  of  the  electron  density 
and  other  quantities.  However,  for  simulations  of  liquids  and  solids,  GGA  XC 
functionals  are  used  in  the  vast  majority  of  DFT  calculations. 

As  noted  above,  DFT  is  currently  the  most  popular  method  for  routine 
quantum  mechanical  calculations.  This  popularity,  which  surged  in  the  early 
1990s,  is  due  to  the  development  of  accurate  XC  functionals.  These  functionals 
offer  accuracy  which  rivals  that  of  high-level  wave  function-based  methods. 
In  terms  of  cost,  DFT  calculations  can  be  performed  on  systems  containing 
up  to  a  few  thousand  atoms  with  modern  computers.  In  terms  of  accuracy, 
studies  [23,24]  show  that  LDA,  BP86  [25,26]  (GGA),  PBE  [27]  (GGA)  and 
B3LYP  [28, 29]  functionals  yield  mean  absolute  errors  of  36.4,  10.3,  8.6  and 
2.2  kca!  mol”1,  respectively,  for  the  atomization  energies  of  the  molecules  in 
the  G-2  set,  [30]  which  is  a  standard  benchmark  database  for  the  assessment 
of  the  accuracy  of  theoretical  methods.  Calculations  of  reaction  barriers  for 
a  set  of  76  reactions  showed  that  LDA,  BP86  [25];  (GGA),  PBE  [27]  (GGA) 
and  B3LYP  [26,29]  (hybrid)  DFT  functionals  yielded  errors  of  14.9,  8.8,  8.7 
and  4.3kcal  mol-1  with  respect  to  the  results  of  high-level  calculations  [31]. 

2.1.3  Plane  Wave  Basis  Sets 

In  most  practical  DFT  calculations  the  Kohn-Sham  single  electron  orbitals, 
are  expressed  as  a  linear  combination  of  basis  functions, 

K 

Mr)  =  ^2c»iXv(r),  (2-8) 

I/=l 

where  v  is  the  index  of  the  basis  function,  K  is  the  total  number  of  basis  func¬ 
tions  and  cVi  is  the  mixing  coefficient  of  basis  function  which  determines 
the  contribution  of  basis  function  Xu  to  orbital  ^7;.  When  optimizing  the  wave 
function,  the  coefficients  cvi  are  variables  that  are  altered  to  yield  the  ‘best5 
wave  function  for  the  system,  while  the  basis  functions  themselves  are  unal¬ 
tered.  The  details  of  the  procedure  used  to  optimize  these  coefficients  will  not 
be  discussed  here;  however,  the  interested  reader  is  directed  to  Chap.  3  of  [32] 
or  Chap.  7  of  [11].  Although  any  type  of  mathematical  function  can  be  used 
as  a  basis  function,  two  types  of  basis  functions  have  achieved  great  popular¬ 
ity:  localized,  atom-centered  functions  and  delocalized  plane  waves.  Localized, 
atom-centered  basis  functions  are  based  on  the  notion  that  a  molecular  system 
is  built  up  from  atoms,  and  hence  the  total  electronic  wave  function  can  be  rep¬ 
resented  as  a  combination  of  atomic  orbital-like  functions,  which  are  centered 
on  specific  nuclei.  Atom-centered  basis  functions  have  achieved  widespread  use 
in  the  quantum  chemistry  community,  where  so-called  Gaussian  functions  [33] 
are  the  most  common  type  of  basis  function.  For  the  work  presented  in  this 
chapter,  localized  basis  functions  are  used  sparingly,  and  will  not  be  discuss 
in  further  detail. 
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An  alternative  to  atom-centered,  localized  basis  functions  involves  the  use 
of  delocalized  functions,  which  are  not  associated  with  specific  nuclei.  The  use 
of  such  functions  is  inspired  by  the  notion  that  electrons  in  periodic  systems 
are  only  slightly  perturbed  by  the  presence  of  nuclei,  and  hence  the  electronic 
wave  functions  should  most  closely  resemble  those  of  free  electrons.  This  de¬ 
scription  is  particularly  suited  to  periodic,  conducting  systems  such  as  metals. 
With  such  systems  in  mind,  one  defines  a  periodic  simulation  cell  (vide  infra), 
e.g.,  a  unit  cell  for  a  crystalline  material,  and  calculates  the  electronic  wave 
function  within  this  cell.  The  natural  basis  set  for  such  a  calculation  is  a  set 
of  plane  waves  [34]: 

Xpw(r)  =  T=eiar,  (2.9) 

where  ft  is  the  volume  of  the  simulation  cell  and  G  is  given  by  the  reciprocal 
lattice  vectors  of  the  cell.  A  plane  wave  basis  set  includes  all  plane  waves  up 
to  some  maximum  value  of  G. 

Plane  wave  basis  sets  have  several  advantages  over  atom-centered  basis 
functions.  First  of  all,  a  plane  wave  basis  set  allows  one  to  move  between  real 
and  reciprocal  spaces  through  the  use  of  Fast  Fourier  Transforms.  This  allows 
one  to  calculate  quantities  in  the  most  convenient  space.  For  example,  deriva¬ 
tives  in  real  space  are  simply  multiplications  in  reciprocal  space,  and  hence  it 
is  convenient  to  perform  such  calculations  in  reciprocal  space.  The  fact  that 
plane  waves  are  not  associated  with  specific  nuclei  also  simplifies  certain  types 
of  calculations.  Along  the  same  lines,  plane  waves  treat  all  regions  of  space 
equally,  and  hence  a  plane  wave  basis  set  is  complete,  which  eliminates  basis 
set  superposition  errors  that  occur  when  using  atom-centered  basis  functions. 
Finally,  plane  waves  are  suited  to  calculations  of  periodic  systems,  which  will 
be  discussed  ahead. 

A  significant  drawback  of  plane  wave  basis  sets  is  that  a  huge  number  of 
plane  waves  must  be  used  to  accurately  describe  the  localized  wavefunctions  in 
regions  around  the  nuclei.  Specifically,  the  wavefunction  and  electron  density 
vary  rapidly  as  one  approaches  a  nucleus.  Accurately  describing  these  rapid 
changes  requires  using  a  plane  wave  expansion  up  to  very  high  values  of  G, 
and  hence  a  large  number  of  plane  wave  basis  functions  must  be  considered 
in  the  calculation.  In  practice,  tens  or  hundreds  of  thousands  of  plane  wave 
basis  functions  must  be  used  to  achieve  the  same  level  of  accuracy  that  can 
be  attained  with  a  few  hundred  Gaussian  basis  functions. 

The  accuracy  and  size  of  the  plane  wave  basis  set  is  determined  by  the 
highest  value  of  G  considered  in  the  plane  wave  expansion.  This  is  determined 
by  taking  advantage  of  the  relationship  between  G  and  the  kinetic  energy  of 
the  plane  wave: 

Ekin{  G)  =  i|G|2,  (2.10) 

which  allows  one  to  define  a  parameter  Ecut  that  corresponds  to  the  kinetic 
energy  of  the  plane  wave  with  the  highest  value  of  G  in  the  basis  set.  The 
parameter  £cu t  is  termed  the  kinetic  energy  cutoff.  Assigning  a  value  for 
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jBcut  indicates  that  all  plane  waves  where  ||G|2  <  Ecut  are  included  in  the 
plane  wave  basis  set.  Without  going  into  details,  this  leads  to  the  following 
expression  for  the  number  of  plane  waves  in  the  basis  set: 

AW  *  ^lElH  (2.11) 

Equation  (2.11)  indicates  that  the  number  of  plane  waves  considered  in  the 
calculation  increases  with  both  the  volume  of  the  simulation  cell,  fi,  and  ECXii . 

Thus  in  summary,  plane  waves  are  an  alternative  to  atom-centered,  local¬ 
ized  basis  sets.  Plane  waves  are  particularly  suited  to  periodic  systems  and 
their  use  simplifies  certain  aspects  of  quantum  chemical  calculations  through 
the  flexibility  offered  by  a  dual  space  representation  and  the  completeness 
of  the  basis  set.  The  major  drawback  of  plane  waves  is  that  a  huge  number  of 
basis  functions  must  be  used  to  attain  a  level  of  accuracy  that  is  comparable 
to  that  achieved  with  even  a  small  number  of  localized  basis  functions. 

A  final  point  that  is  closely  related  to  basis  functions  involves  the  use  of 
pseudopotentials  (also  known  as  effective  core  potentials)  to  replace  the  core 
orbitals  and  to  smooth  the  valence  orbitals  near  the  nuclei.  For  example,  in 
carbon,  which  has  a  ls22s22p2  electronic  configuration,  the  2s  and  2p  electrons 
would  be  treated  with  a  basis  set,  while  the  Is  electrons  would  be  represented 
by  a  pseudopotential.  The  motivation  for  eliminating  the  explicit  treatment 
of  core  electrons  stems  from  the  notion  that  these  electrons  do  not  contribute 
significantly  to  chemical  bonding,  and  hence  the  core  states  are  altered  min¬ 
imally  during  chemical  reactions.  As  such,  the  core  states  do  not  need  to  be 
treated  explicitly  in  most  quantum  mechanical  calculations.  This  dramatically 
decreases  the  number  of  basis  functions  needed  for  a  particular  calculation, 
which  decreases  computational  cost.  An  additional  advantage  of  pseudopoten¬ 
tials  is  that  relativistic  effects  associated  with  core  electrons  in  heavy  atoms 
can  easily  be  incorporated  into  the  calculation  without  additional  computa¬ 
tional  effort.  In  practice,  all  plane  wave  calculations  use  pseudopotentials  to 
represent  the  core  states  of  all  atoms. 

2.1.4  Periodic  Boundary  Conditions 

It  is  often  of  interest  to  perform  simulations  of  bulk  liquids  or  bulk  solids. 
By  definition,  such  systems  extend  across  length  scales  that  are  inaccessible 
through  quantum  mechanical  (or  even  classical)  simulation.  One  means  of 
simulating  these  systems  is  to  perform  a  calculation  on  a  small  system  that  is 
representative  of  the  bulk  material  and  repeat  this  small  system  infinitely  in 
all  three  spatial  directions.  A  two-dimensional  example  is  shown  in  Fig.  2.3, 
where  the  central  cell  is  repeated  in  the  plane  of  the  page.  Representing  a 
system  in  this  manner  is  termed  using  periodic  boundary  conditions  (PBCs). 
In  what  follows,  a  few  key  points  regarding  the  use  of  PBCs  in  quantum 
mechanical  calculations  are  qualitatively  described. 
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Fig.  2.3.  Periodic  boundary  conditions  in  two  dimensions.  The  central  cell  contains 
five  water  molecules  and  is  repeated  exactly  in  both  spatial  directions  that  define 
the  plane  of  the  page.  Nine  cells  are  shown  explicitly;  however,  the  system  should 
be  repeated  infinitely  in  all  directions 


As  mentioned  above,  the  use  of  PBCs  involves  defining  a  simulation  cell 
and  repeating  this  cell  periodically  to  form  an  infinite  crystal.  In  this  scheme, 
quantities  such  as  the  wave  function  are  periodic  and  must  be  continuous 
across  the  walls  of  neighboring  simulation  cells.  This  introduces  constraints 
into  the  calculations  which  ensure  that  the  central  simulation  cell  inter¬ 
acts  with  its  periodic  images,  although  the  calculation  is  only  performed  on 
the  central  simulation  cell.  Long  range  electrostatic  interactions  between  the 
periodic  images,  can  be  accounted  for  by  using  techniques  such  as  Ewald 
summations  [35].  This  property  of  PBCs  proves  useful  when  simulating  bulk 
materials,  where  any  given  atom  should  interact  with  those  throughout  the 
material.  For  example,  in  a  crystal,  which  is  composed  of  periodically  repeated 
unit  cells,  the  atoms  in  a  given  unit  cell  should  interact  with  those  in  the  other 
unit  cells.  In  the  case  of  non-crystalline  systems,  e.g.,  liquids  or  amorphous 
solids,  it  is  desirable  to  study  bulk  systems;  however,  such  systems  should  not 
be  periodic.  PBCs  can  still  be  used  to  approximate  these  systems  by  defining 
a  large  “central”  simulation  cell,  which  minimizes  the  effect  of  externally  im¬ 
posed  periodicity  on  the  behavior  of  the  system.  It  should  be  noted  that 
imposing  PBC  implies  some  constraints  on  the  long-range  structure  and  cor¬ 
relations  in  bulk  systems.  These  must  be  carefully  considered  when  simulating 
long-range  properties  such  as  low-frequency  phonon  modes. 

2.1.5  Molecular  Dynamics 

In  the  preceding  section,  we  have  outlined  how  first  principles  density  fuin> 
tional  theory  calculations  can  be  used  to  calculate  the  potential  energy  HMP* 
face.  Exploring  the  PES  is  important  because  it  provides  insight  into  OhOfiiical 
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processes  that  may  occur.  Conventional  studies  of  materials  and  molecules 
that  apply  first  principles  quantum  mechanical  methods  often  do  so  in  a  sta¬ 
tic  sense.  That  is,  the  behavior  on  the  PES  is  considered  in  the  zero  Kelvin 
limit  and  then  extrapolated  to  finite  temperatures.  Furthermore,  static  cal¬ 
culations  focus  on  particular  regions  of  the  PES  corresponding  to  relevant 
species  such  as  the  presumed  reactants,  products  and  transition  state  associ¬ 
ated  with  a  given  processes.  ‘Static’  calculations  are  extremely  powerful  and 
widely  used.  In  this  chapter,  more  novel  ab  initio  molecular  dynamics  (MD) 
simulations  are  also  presented.  These  ab  initio  MD  simulations  allow  one  to 
observe  the  evolution  of  a  system  in  space  and  time  at  the  molecular  level 
subject  to  a  given  set  of  macroscopic  external  conditions,  e.g.,  temperature 
and  pressure.  It  is  noted  that  MD  simulations  provide  a  great  deal  of  quali¬ 
tative  insight  into  chemical  processes,  yet  can  also  be  applied  quantitatively. 
In  what  follows,  we  discuss  the  MD  methodology  and  its  use  in  the  study  of 
chemical  reactions. 

Most  MD  calculations  treat  the  nuclei  comprising  a  chemical  system  as 
classical  particles  whose  trajectories  can  be  determined  with  Newtonian  equa¬ 
tions  of  motion: 

d2Ri 

F  t  =  rmaa  =  m*  — 2~  (2.: 12) 

where  the  force  on  atom  i,  Ft,  is  obtained  by  differentiating  the  potential 
energy  with  respect  to  the  nuclear  coordinates: 


dE(  R) 
dR  i 


(2.13) 


Unfortunately,  it  is  not  possible  to  analytically  integrate  (2.12)  for  a  many- 
body  system,  and  one  must  numerically  integrate  the  system  of  coupled  differ¬ 
ential  equations  representing  the  trajectory  over  a  series  of  small  time  steps. 
The  general  procedure  is  as  follows.  At  a  certain  time,  t ,  the  nuclei  are  de¬ 
scribed  by  a  set  of  positions  and  momenta,  R*(0)  and  m7;V*(0).  The  potential 
energy  for  the  system  is  calculated  and  used  to  determine  the  forces  on  the 
nuclei,  which  are  then  used  to  update  the  nuclear  positions  and  momenta  to 
their  values  at  a  later  time,  t  +  At.  The  process  is  then  repeated  until  a  suffi¬ 
ciently  long  total  period  of  time  is  simulated.  The  time  intervals,  At ,  used  in 
this  process  must  be  short  enough  to  capture  the  highest  frequency  motions 
in  the  system,  which  typically  correspond  to  bond  vibrations.  As  such,  these 
intervals  are  on  the  order  of  10-]5s,  and  hence  ~106  force  calculations  must 
be  performed  to  simulate  one  nanosecond  of  time.  Despite  the  considerable 
computational  expense  associated  with  this  procedure,  MD  simulation  has 
been  an  exceedingly  powerful  and  popular  tool  for  over  50  years  [36-38]. 

An  important  feature  of  MD  simulations  is  the  ability  to  study  reactions 
under  a  given  set  of  conditions.  For  instance,  by  properly  scaling  the  nuclear 
kinetic  energies  it  is  possible  perform  simulations  at  a  given  temperature. 
Means  of  applying  pressure  to  the  system  also  exist.  In  a  constant  volume  sim¬ 
ulation  with  PBC,  a  simulation  with  high  molecular  densities  will  effectively 
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lead  to  high  simulation  pressures.  Alternatively,  pressure  and  temperature  can 
be  added  as  independent  dynamic  variables  in  the  equations  of  motion  that  arc 
coupled  to  the  positions  and  momenta  of  the  molecules  in  the  simulation  [35] . 
These  capabilities  are  quite  important  when  studying  the  influence  of  external 
conditions  on  a  reaction.  By  comparison,  temperature  and  pressure  can  only 
be  incorporated  into  static  calculations  in  an  ad  hoc  manner. 

2.1.6  Ab  initio  Molecular  Dynamics 

In  the  preceding  subsection,  the  classical  MD  methodology  was  briefly  in¬ 
troduced  in  a  generic  sense,  with  the  PES  considered  merely  as  a  means  of 
obtaining  the  forces  on  the  nuclei.  For  the  vast  majority  of  MD  simulations, 
the  PES  is  evaluated  with  classical  potentials,  which  are  parameterized  func¬ 
tions  that  subdivide  the  potential  energy  into  distinct  contributions  associated 
with  bond  lengths,  angles,  dihedrals,  electrostatics,  etc.  With  the  exception  of 
a  few  specific  cases,  classical  potentials  cannot  accurately  describe  changes  in 
bonding,  and  hence  are  not  suitable  for  use  in  the  simulation  of  chemical  reac¬ 
tions,  i.e.,  bond  formation  and  dissociation.  In  order  to  study  such  processes, 
it  is  necessary  to  explicitly  consider  the  electronic  structure  of  the  system. 
A  popular  way  to  do  this  is  ab  initio  molecular  dynamics  (AIMD)  [39]. 

In  the  AIMD  approach,  the  potential  energy  in  (2.13)  is  derived  from  a 
QM  calculation;  typically  at  the  density  functional  level  of  theory.  The  ex¬ 
plicit  consideration  of  the  electronic  structure  of  the  system  in  this  manner 
leads  to  a  powerful  simulation  method  that  is  capable  of  describing  chemical 
reactions.  In  fact,  because  the  system  moves  on  the  electronic  PES  accord¬ 
ing  to  well-defined  equations  of  motion,  i.e.,  Newton’s  equations,  one  can,  in 
principle,  simply  initiate  an  AIMD  simulation,  follow  the  motion  of  the  sys¬ 
tem  and  observe  chemical  reactions  without  having  any  preconceived  notions 
of  what  the  potential  energy  surface  describing  those  reactions  may  be.  The 
ability  of  AIMD  simulations  to  investigate  known  reactions  and  identify  new 
reactions  has  led  to  a  dramatic  increase  in  the  popularity  of  this  method  in 
recent  years. 

AIMD  simulation  is  not  without  its  drawbacks,  as  these  calculations  suffer 
from  many  limitations  due  to  the  computational  expense  of  the  first  principles 
QM  calculations.  If  the  potential  energy  is  derived  from  a  DFT  calculation, 
it  is  typically  only  possible  to  perform  simulations  on  systems  composed  of 
a  few  hundred  atoms,  and  only  sub- nanosecond  time  scales  may  be  consid¬ 
ered.  Despite  these  limitations,  AIMD  simulation  has  recently  experienced  a. 
tremendous  increase  in  popularity  and  is  commonly  applied  to  a  wide  variety 
of  systems  [4, 5,40-42]. 

Two  popular  types  of  AIMD  simulations  exist.  The  first  approach  is  called 
Born-Oppenheimer  AIMD  (BO- AIMD)  and  involves  optimizing  the  electronic 
structure  of  the  system  at  each  time  step  of  the  simulation.  That  is,  a  full 
wave  function  optimization  is  performed  at  each  time  step,  and  the  resulting 
potential  energy  is  used  to  determine  the  forces  on  the  nuclei.  BO-AIMI)  is 
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an  intuitive  approach  to  AIMD  simulation,  yet  is  computationally  expensive 
because  of  the  tremendous  number  of  wave  function  optimizations  that  must 
be  performed. 

The  second  main  type  of  AIMD  simulation  is  known  as  Car-Parrinello 
AIMD  (CP-AIMD)  [3]  and  involves  treating  the  electronic  states  as  effective 
dynamic  variables  which  are  propagated  according  to  the  classical  Newtonian 
equations  of  motion.  In  a  CP-AIMD  calculation,  the  first-principles  method 
used  is  almost  exclusively  DFT,  where  the  Kohn-Sham  orbitals  are  optimized 
at  the  very  first  step  of  the  simulation  only.  Then,  they  are  assigned  a  fictitious 
mass,  /i,  and  propagated  in  time  according  to  classical  equations  of  motion 
through  the  following  extended  Lagrangian: 


-t  N atom  a  1  ^lorbttais 

CCP  =  \  E  +  \  E  Mi  (xlH\ih)  -  Edft\p) 


7  =  1 
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orthonormality  constraints 


(2.14) 


The  first  term  of  (2.14)  is  the  classical  kinetic  energy  of  nuclei  with  masses  Mj, 
the  second  term  is  the  fictitious  ‘classical’  kinetic  energy  of  orbitals  due  to  the 
electronic  degrees  of  freedom  with  inertial  parameters  (fictitious  masses)  pi9 
the  third  term  is  the  electronic  energy  as  described  in  (2.7),  and  the  last  term 
is  the  constraints  to  maintain  orthogonality  of  orbitals.  Upon  substitution 
of  (2.14)  into  Lagrange’s  equation  of  motion  one  obtains  the  Car-Parrinello 
equation  of  motion  for  both  the  nuclei  and  the  KS  orbitals: 


= 


MjRj  : 

dEtyuIh) 
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_  dEj'ipj,  flj) 

dlh 

+  orthogonolity  constraints. 


(2.15) 

(2.16) 


In  this  approach,  the  Kohn-Sham  orbitals  are  propagated  as  classical  degrees 
of  freedom  at  the  same  time  the  nuclear  motion  is  treated.  Therefore,  there 
is  no  need  to  optimize  the  electronic  structure  at  each  time  step  as  in  the 
BO-AIMD  approach.  It  is  important  to  note  that  although  the  orbitals  are 
treated  as  classical  degrees  of  freedom,  the  electronic  structure  is  still  deter¬ 
mined  quantum  mechanically.  It  is  also  noted  that  the  mass,  //,  is  a  fictitious 
quantity  assigned  to  the  orbitals  that  should  not  be  confused  with  the  mass  of 
the  electron.  Moreover,  in  the  limit  that  p  approaches  zero,  the  Car-Parrinello 
equation  converges  to  the  KS  equation  so  that  E  is  minimized  and  the  Born- 
Oppenheimer  limit  is  recovered. 

The  key  -to  successfully  performing  a  CP-AIMD  simulation  lies  in  prevent¬ 
ing  the  transfer  of  energy  between  the  nuclear  and  electronic  subsystems.  This 
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is  achieved  by  altering  the  value  of  /x  such  that  the  power  spectra  of  these  two 
subsystems  do  not  overlap,  i.e.,  \x  should  be  kept  small  to  ensure  that  the 
orbitals  oscillate  at  frequencies  much  higher  than  those  of  the  nuclei.  Further¬ 
more,  using  a  small  value  for  fj,  allows  the  orbitals  to  adjust  quickly  to  changes 
in  the  nuclear  configuration,  which  keeps  the  system  close  to  the  ground  state 
PES.  In  a  properly  performed  CP-AIMD  simulation  the  electronic  structure 
will  remain  close  to  the  true  PES,  and  hence  the  dynamics  observed  in  a 
CP-AIMD  simulation  closely  follow  those  which  would  be  determined  in  a 
BO-AIMD  calculation. 

As  indicated  above,  the  main  advantage  of  the  CP-AIMD  approach  is 
that  the  electronic  structure  is  only  optimized  once  during  the  course  of  a 
simulation.  This  can  significantly  decrease  computational  cost  compared  with 
a  BO-AIMD  simulation,  where  the  electronic  structure  is  optimized  at  each 
time  step.  However,  properly  integrating  the  orbital  dynamics,  which  oscillate 
at  high  frequencies,  requires  the  use  of  times  steps,  At,  that  are  typically 
an  order  of  magnitude  smaller  than  those  used  in  a  BO-AIMD  simulation. 
This  decreases  the  computational  advantage  of  CP-AIMD  simulations  over 
BO-AIMD  simulations.  Nonetheless,  the  former  is  faster  in  most  cases  and  has 
become  the  most  popular  AIMD  simulation  technique  for  studying  chemical 
systems. 

2.2  Collision  Dissociation  of  Nitromethane 

Although  numerous  theoretical  reports  focusing  on  nitromethane  (NM)  have 
appeared  in  the  literature  [43-50]  the  decomposition  mechanism  of  liquid  NM 
under  shock  conditions  remains  elusive.  The  C — N  bond  is  the  weakest  bond  in 
NM  and  therefore  unimolecular  C — N  bond  scission  has  often  been  implicated 
as  the  dominant  decomposition  pathway.  However,  a  number  of  high  pressure 
studies  do  not  support  this  mechanism,  [51,52]  vide  infra.  To  gain  some  insight 
into  these  decomposition  pathways  bimolecular  collision  MD  simulations  have 
been  performed  by  a  number  of  researchers  [53,54].  A  summary  [54]  of  bimole¬ 
cular  collisions  of  NM  in  a  variety  of  orientations  and  a  large  range  of  incident 
collision  velocities  is  shown  in  Table  2.1.  For  the  offset  anti-parallel  molecular 
orientation,  the  critical  velocity  for  successful  dissociation  was  found  to  be 
7.0km  s"\  higher  than  the  average  atomic  velocities  expected  at  the  shock 
front  of  a  detonation  [55].  Although  these  bimolecular  collision  simulations 
provide  a  qualitative  glimpse  into  the  shock-induced  dissociation  process,  they 
cannot  account  for  the  effects  that  neighboring  molecules  may  have  on  the  dm° 
sociation  mechanism.  This  issue  is  addressed  in  first  principles  Car  Parriuollo 
molecular  dynamics  (CPMD)  simulations  of  multimolecular  collisions  ol  NM. 

2.2.1  Impact  of  a  Single  Molecule  on  Multiple  Molecules 

The  previous  bimolecular  collision  simulations  are  limited  to  collision-induced 
reactions  without  neighboring  confinement  and  collision  sequences  beyond 
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Table  2.1.  Threshold  collision  velocities  (km  s"*1)  determined  from  multimolecular 
and  bimolecular  collision  simulations 

Orientation  Threshold  collision  velocity 


Multimolecular  Bimolecular 


H\  /? 

HH  O 

o^-so 

N 

/’SgH 

H  H 

Perpendicular 

12.0 

8.0 

HV 

“Y  ® 
0^0 

O^o 

N 

/SjlH 

H  H 

Anti-parallel 

11.0 

10.5 

h\_/9  . 

HH  O 

HH  0 

Head-to-tail 

11.0 

8.5 

ay 

0^0  ^ 

>4 

IX  O 

Offset  Anti-parallel 

12.0 

7.0 

% 

the  initial  impact  are  meaningless  due  to  only  the  two  isolated  molecules  in¬ 
volved.  Multimolecular  collision  simulations  offer  a  substantial  improvement 
since  they  allow  for  possible  collision-induced  reactions  involving  neighboring 
molecules  beyond  the  initial  impact.  The  effects  of  multiple  collisions  on  the 
decomposition  mechanism  are  important  in  that  they  determine  whether  the 
molecular  decomposition  occurs  at  the  initial  shock  wave  front  or  after  further 
thermalization  of  the  shock  wa  ve  energy.  This  issue  can  have  important  impli¬ 
cations  on  the  behavior  of  the  material  as  a  high  energy  density  material.  Due 
to  the  very  short  time  scales  involved  in  passing  of  the  initial  shock  wave  and 
thermalization  stages,  information  regarding  the  decomposition  mechanism  at 
these  time  scales  can  be  experimentally  difficult  to  obtain.  This  makes  the  use 
of  simulation  methods  very  useful  for  the  study  of  this  phenomenon. 

The  first  study  of  multimolecular  collisions  involves  collisions  by  the  im¬ 
pact  of  a  single  molecule  on  multiple  molecules.  The  simulations  are  conducted 
under  the  same  four  collision  orientations  of  the  initial  colliding  NM  pair  of 
molecules  as  the  bimolecular  collisions  reported  in  Table  2.1.  A  total  of  13  NM 
molecules  are  placed  in  an  orthorhombic  simulation  cell  of  19.0  x  14.2  x  14.2  A3 
such  that  one  molecule  is  placed  at  one  end  of  the  simulation  cell  and  the 
remaining  12  molecules  are  placed  at  the  opposite  end.  The  simulations  are 
initiated  by  propelling  the  separated  molecule  into  a  stationary  molecule, 
embedded  in  the  cluster  of  NM  molecules  as  illustrated  in  Fig.  2.4.  The  cluster 
of  12  molecules  is  randomly  oriented  and  possesses  a  density  of  1.14 g  cm*3, 
corresponding  to  the  density  of  liquid  NM  at  ambient  pressure  and  temper¬ 
ature.  The  incident  collision  velocities  range  from  6  to  12  km  for  various 
orientations  of  the  two  colliding  molecules.  Each  multimolecular  collision  sim¬ 
ulation  is  run  for  3,500  time  steps  or  0.507  ps. 
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Fig.  2.4.  Schematic  of  the  single-molecule  on  multi- molecule  nitromethane  collision 
simulation  with  the  head-to-tail  orientation  for  the  initially  colliding  pair 


As  evident  from  Table  2.1,  the  threshold  collision  velocities  determined 
from  the  multimolecular  collision  simulations  are  significantly  higher  than 
those  found  previously  from  bimolecular  collision  simulations  for  each  mole¬ 
cular  orientation  investigated.  In  all  cases  of  multimolecular  collisions,  the 
sub-critical  velocity  collisions  induce  initial  C  N  bond  scission  of  the  in¬ 
coming  and  stationary  molecules.  However,  the  neighboring  molecules  act  as 
a  trap  to  confine  the  CH3  and  NO2  fragments  produced,  thereby  enabling 
them  to  recombine  to  form  intact  NM  molecules.  The  C  N  bond  then  alter¬ 
nately  breaks  and  reforms  as  the  simulation  continues  until  finally  it  remains 
intact  when  there  is  insufficient  energy  to  break  it.  At  the  threshold  colli¬ 
sion  velocities,  the  recoiling  fragments  produced  during  the  initial  collision 
possess  enough  translational  energy  to  overcome  the  confinement  forces  of 
the  neighboring  molecules  leading  to  permanent  cleavage  of  the  C  N  bond 
to  yield  CH3  and  NO2.  Collisions  of  more  than  two  molecules  are  observed, 
for  instance,  for  the  perpendicular  molecular  orientation  at  12  km  s-1,  where 
collisions  occur  between  four  NM  molecules  including  the  incoming,  the  sta¬ 
tionary  and  two  neighboring  molecules.  This  results  in  permanent  C — N  bond 
cleavage  of  the  stationary  molecule,  while  the  other  molecules  involved  recom¬ 
bine  to  form  intact  NM  molecules.  Permanent  C — N  bond  scission  is  observed 
in  both  the  incoming  and  stationary  NM  molecules  in  the  offset  anti-parallel 
orientation  at  its  threshold  velocity  of  12  km  s-1,  while  the  anti-parallel  and 
head-to-tail  orientations  display  permanent  C — N  bond  cleavage  of  the  in¬ 
coming  NM  molecule  at  the  threshold  velocity  of  11  km  s-1.  These  velocities 
are  not  available  under  conventional  detonation  conditions  of  nitromethane 
with  a  shock  velocity  of  6-7  km  s-1  where  the  molecular  collision  velocities 
are  even  less  than  this  value. 
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While  the  threshold  velocity  of  multimolecular  collisions  is  higher  than 
that  of  bimolecular  collisions  in  all  the  collision  orientations,  the  predominant 
mechanism  of  decomposition  remains  the  C — N  bond  cleavage  as  in  the  bi¬ 
molecular  collisions.  A  different  fragmentation  mechanism,  initial  C  H  bond 
scission  of  the  stationary  molecule  is  observed  in  one  of  the  high  velocity 
multimolecular  collision  simulations  (12  km  s-1  in  a  head-to-tail  orientation). 
The  initial  C — H  bond  cleavage  is  followed  immediately  by  the  migration  of 
this  free  H  atom  to  one  of  the  O  atoms  of  the  incoming  NM  molecule  which 
induces  N  O  bond  cleavage  of  this  NM  molecule  yielding  OH  and  CH3NO. 
The  N — O  bond  of  the  stationary  molecule  is  then  cleaved  and  the  free  O 
atom  migrated  to  the  CH2  group  of  the  stationary  molecule,  thereby  inducing 
C  N  bond  cleavage  to  yield  CH20  and  NO  fragments.  The  final  products 
produced  from  this  cascade  of  reactions  are:  OH,  CH20,  NO,  and  CH3NO. 

2.2.2  Impact  of  Multiple  Molecules  on  Multiple  Molecules 

The  multimolecular  collision  simulations  are  further  extended  to  the  collisions 
between  two  groups  of  molecules  as  depicted  in  Fig.  2.5.  The  chemical  system 
of  32  NM  molecules  is  contained  in  a  cubic  simulation  cell  14.2  A  on  each  side. 
The  density  of  the  32  molecules  in  the  cube  corresponds  to  1.14  g  cm-  *  for 
liquid  NM.  Periodic  boundary  conditions  are  employed  to  better  represent  the 
bulk  material  properties  and  effects  of  neighboring  molecules.  In  an  identical 
setup  without  using  periodic  boundary  conditions,  molecules  leaving  the  cell 
will  effectively  be  in  a  vacuum  and  eventually  disperse  as  a  result  of  Brownian 
motion.  Initial  positions  for  the  simulation  were  first  obtained  from  a  classical 
MD  simulation.  The  initial  orientations  and  velocities  of  each  of  the  NM  mole¬ 
cules  in  the  simulation  cell  are  derived  from  an  AIMD  of  a  sample  that  was 
equilibrated  for  approximately  4.36  ps  at  298.15  K  with  a  Nose- Hoover  ther¬ 
mostat.  Eight  NM  molecules,  occupying  the  most  right-hand  quarter  of  the 
simulation  cell,  are  translated  linearly  forward  to  impact  into  the  other  24  NM 
molecules  occupying  the  rest  of  the  cell  as  shown  in  Fig.  2.5.  In  each  simu¬ 
lation,  all  eight  NM  molecules  are  assigned  a  uniform  initial  velocity  ranging 
from  2  to  10  km  s-1. 

The  simulations  are  further  characterized  as  static  or  dynamic  simulations. 
The  static  simulations  are  similar  to  the  previous  simulations  of  bimolecular 
collisions  and  single-to-multiple  molecular  collisions,  in  which  the  kinetic  en¬ 
ergy  and  temperature  of  the  molecules  are  initially  set  to  zero.  The  eight  mole¬ 
cules  occupying  the  right-hand  quarter  of  the  cube  are  then  translated  forward 
with  their  assigned  velocity.  In  the  dynamic  simulations,  all  32  molecules  ini¬ 
tially  possess  equilibrium  velocity  and  the  total  kinetic  energy  corresponds  to 
the  kinetic  energy  of  the  equilibrated  sample  of  liquid  NM  at  298.15  K  and 
1  atm.  The  eight  molecules  occupying  the  right-hand  quarter  of  the  simulation 
cube  are  translated  with  an  assigned  velocity  in  addition  to  their  equilibrium 
velocity.  Thus,  the  dynamic  simulations  are  considered  more  realistic  than 
the  static  simulations.  All  simulations  were  performed  for  5,000  time  steps  or 
0.725  ps. 
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Fig.  2.5.  Schematic  of  the  multiple-molecule  on  multiple-molecule  collision 
simulations 


Figure  2.6  shows  the  dynamic  simulation  results  for  average  C  N  bond 
lengths  of  all  32  NM  molecules  as  a  function  of  time  at  collision  velocities 
ranging  from  2  to  10  km  s_1.  At  lower  velocities,  the  bond  lengths  oscillate 
around  the  initial  bond  lengths.  The  initial  bond  length  at  25° C  is  slightly 
higher  than  its  ideal  bond  length,  which  is  defined  for  0  K  in  the  gas  phase  and 
imparts  the  lowest  total  energy  to  the  molecule.  The  ideal  bond  lengths  for 
the  C  N,  N  O  and  C  H  are  1.46,  1.22  and  1.10  A,  derived  from  geometry 
optimizations  of  the  isolated  nitromethane  molecule  at  the  PBE/6-3lG(d,p) 
level.  A  bond  is  considered  dissociated  if  the  distance  between  the  atoms  of 
the  bond  is  greater  than  1.5  times  the  ideal  bond  length  which  corresponds 
to  2.19  A  for  the  C — N  bond.  No  bond  breakage  was  observed  at  velocities 
of  2-6  km  s-1.  However,  for  a  collision  velocity  of  8  km  s-1,  there  is  a  slight 
increase  in  the  average  C — N  bond  distance  in  the  simulations  approximately 
500  time  steps  into  the  simulation,  which  further  increases  at  around  1,250 
time  steps.  Qualitative  examination  of  the  simulation  reveals  that  the  C — N 
bond  of  a  single  nitromethane  molecule  dissociated  at  approximately  1,250 
time  steps.  However,  the  neighboring  molecules  confine  the  CH3  and  NO2 
fragments  produced,  thereby  enabling  them  to  recombine  to  form  an  intact 
NM  molecule  at  nearly  1,600  time  steps.  The  C — N  bond  then  alternately 
breaks  and  reforms  as  the  simulation  continues  until  finally  it  remains  intact 
when  there  is  insufficient  kinetic  energy  to  overcome  the  bond  strength  and  the 
confinement  forces  of  neighboring  molecules.  At  8  km  s_1,  while  the  average 
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Fig.  2.6.  Average  C — N  bond  length  as  a  function  of  time  at  2-10  km  s-1  velocity 
of  impact  of  multiple  molecules  on  multiple  molecules.  The  initial  temperature  and 
pressure  of  the  molecules  in  the  simulation  are  25° C  and  1  atm 


N — O  and  C  H  bond  lengths  oscillate  around  their  initial  bond  lengths  at 
enlarged  fluctuation  amplitude,  no  N  O  and  C  H  bonds  are  broken. 

As  the  collision  velocity  is  increased  to  10  km  s'1  the  first  C  N  bond 
breakage  is  observed  at  610  time  steps,  resulting  in  an  increase  in  average 
bond  length  to  1.6  A  (Fig.  2.6).  The  time  delay  from  the  initial  impact  is 
~0.05ps,  which  is  shorter  than  the  time  delay  at  the  collision  velocity  of 
8  km  s"1  (~0.13ps).  The  time  delay  between  the  initial  impact  and  the  C  N 
bond  breakage  and  its  dependence  on  impact  velocity  are  likely  associated 
with  the  conversion  of  translational  energy  to  vibrational  energy  required  for 
cleavage  of  the  bond.  As  time  proceeds,  more  and  more  C  N  bonds  are 
permanently  broken  as  shown  in  Fig.  2.7.  Breakage  of  the  N  O  and  C — H 
bonds  occurs  only  at  10  km  s-1  and  in  a  much  delayed  time.  Nearly  2,600  time 
steps  elapse  till  the  first  N  O  bond  dissociation  and  3,000  time  steps  elapse 
till  the  first  C  H  bond  dissociation.  In  all  cases  of  N  O  bond  dissociation, 
bond  cleavage  occurs  on  an  intact  NM  molecule  and  results  in  migration  of 
the  oxygen  atom  to  a  methyl  fragment  forming  a  formaldehyde  intermediate. 
Also,  the  interaction  of  two  methyl  fragments  repeatedly  forms  and  breaks 
short-lived  ethane.  The  C  H  bond  cleavage  results  in  only  a  small  number  of 
hydrogen  atoms  that  migrate  back  and  forth  to  the  nitro  group  of  neighboring 
nitromethane. 

The  fact  of  the  earlier  C — N  bond  cleavage  suggests  that  the  predomi¬ 
nant  mechanism  of  decomposition  remains  the  C — N  bond  breakage  into  the 
nitro  and  methyl  groups  at  the  threshold  collision  velocity  between  8  and 
10  km  s-1.  As  indicated  in  Table  2.1,  this  threshold  velocity  is  less  than  the 
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Fig.  2.7.  Number  of  C — N  bonds  broken  as  a  function  of  time  at  10  km  s_1  velocity 
of  impact  of  multiple  molecules  on  multiple  molecules.  The  initial  temperature  and 
pressure  of  the  molecules  in  the  simulations  are  25° C  and  1  atm 


values  obtained  from  the  collisions  of  a  single  molecule  on  multiple  molecules. 
This  is  so  because  in  the  collisions  of  multiple  molecules  on  multiple  molecules, 
the  total  impact  kinetic  energy  received  by  the  receptor  molecules  is  larger 
and  the  latent  neighboring  confinement  becomes  weaker  due  to  the  distributed 
impact  kinetic  energy. 

The  static  simulations  provide  an  identical  decomposition  mechanism  and 
threshold  collision  velocity  based  on  the  C  N  bond  cleavage.  The  difference 
between  the  static  and  dynamic  simulations  lies  in  the  time  delay  of  the  de¬ 
composition  and  the  rate  of  bond  length  increase  that  becomes  apparent  at  the 
impact  velocity  of  10  km  s-1  as  shown  in  Fig.  2.6.  The  difference  observed  can 
be  explained  by  the  additional  kinetic  energy  at  room  temperature  and  pres¬ 
sure  available  in  the  dynamic  simulations.  This  additional  energy  manifests 
itself  vibrationally  and  promotes  decomposition  once  initial  bond  breakage 
begins. 

It  is  clear  from  these  simulations  that  a  strong  shock  can  result  in  rapid 
decomposition  of  nitromothane,  in  times  far  short  for  thermal  equilibrium  to 
have  been  established.  The  simulations  of  this  section  indicate  t  hat  the  ini¬ 
tial  velocities  required  for  the  direct  decomposition  of  nitromet  bane  molecules 
is  in  the  range  of  8  10  km  s”1.  Experimentally,  it,  was  found  that  the  aver¬ 
age  speed  of  nitromethane  molecules  at  the  detonation  shock  front  is  about 
2.5km  s"1  [56].  This  shows  that  the  nitromethane'  decomposition  react  ion 
mainly  occurs  in  the  thermalization  stage  behind  the  shock  front  of  a  con¬ 
ventional  detonation.  This  conclusion  also  appears  to  be  continued  by  recent 
experiments  of  the  incident  or  multiple  reflected  shock  compression  initiation 
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of  detonation  in  li<|iii« I  nil  romcthane,  where  the  temperatures  for  incident  and 
multiple  reflected  initiation  upon  shock  compression  arc  equal  to  about  700  K, 
while'  the  multiple  reflected  shock  pressure  reaches  8GPa,  nearly  two  times 
that  of  the  incident  shock  [57,58].  The  8  10km  s-1  decomposition  threshold 
collision  velocity  further  suggests  that  the  collision-induced  detonation  initi¬ 
ation  in  the  shock  front  would  likely  be  dominant,  for  a  detonation  velocity 
beyond  15  20  km  s-1  in  a  molecular  condensed  matter.  This  would  result  in  a 
megabar  detonation  pressure  since  the  detonation  pressure  is  approximately 
proportional  to  the  square  of  the  detonation  velocity.  In  order  to  examine  the 
threshold  pressure  for  mechanical  dissociation  under  minimal  thermal  influ¬ 
ence,  first  principles  molecular  dynamics  simulations  of  static  compression  at 
low  temperatures  will  be  discussed  in  the  next  section. 


2.3  Pressure  Dissociation  of  Nitromethane 

The  reactions  of  high  energy  density  materials  under  low  temperature  high 
pressure  conditions  are  of  interest.  For  this  purpose,  the  high  pressure  decom¬ 
position  of  nitromethane  has  been  studied.  To  study  the  effects  of  isotropic 
pressure  on  the  decomposition  of  NM,  Car  Parrinello  ab  initio  MD  simu¬ 
lations  were  performed.  Periodic  boundary  conditions  were  utilized  with  a 
simple  cubic  simulation  cell  and  32  nitromethane  molecules  such  that  the 
experimentally  measured  bulk  density  of  1.139 gem-3  at  298 K  and  1  atm 
pressure  was  reproduced.  This  simulation  cell  volume  was  then  shrunk  by  a 
factor  of  1.5  3.0.  Well  equilibrated  classical  MD  simulations  were  used  as  the 
starting  point. 

Initial  geometric  configurations  for  the  ab  initio  MD  simulations  were  ex¬ 
tracted  from  well  equilibrated  classical  MD  simulations.  The  cell  volume  was 
compressed  isotropically  by  a  factor  of  up  to  3.  Slow  heat  up  was  performed 
(50 K  — *  150 K  — *  3(H) K).  At  300 K,  equilibration  was  performed  for  25,000 
time  steps  (ea.  3.5  ps)  and  production  runs  were  performed  for  25,000  time 
steps.  The  Ron,  KN().  and  RCH  bond  lengths  at  different  compressions  are 
shown  in  Fig.  2.N. 

The  ab  initio  MD  simulations  of  compressed  liquid  NM  show  no  molecular 
decomposition  at.  compression  factors  of  2.0  and  2.5  at  300  K.  However,  when 
the  liquid  is  compressed  by  a  factor  of  3.0  an  interesting  pressure-induced 
molecular  decomposition  is  observed  at  150  K  after  equilibration  for  about 
3,300  time  steps.  An  intormolecular  H  atom  transfer  event  occurs  between 
two  very  closely  spaced  molecules  aligned  in  an  anti-parallel  alignment,  as 
illustrated  in  the  MD  trajectory  snapshots  shown  in  Fig.  2.9.  As  illustrated  in 
fig.  2.9,  the  CH3  group  of  one  of  the  NM  molecule  re-orients  itself  such  that 
two  of  its  H  atoms  lie  in  the  same  plane  as  the  N02  fragment  and  the  other 
H  atom  is  oriented  nearly  perpendicular  to  this  plane  directed  at  the  ()  atom 
of  the  adjacent  NM  molecule.  This  H  atom  then  transfers  from  the  C  atom  of 
the  first  NM  molecule  to  the  ()  atom  of  the  second  NM  molecule  to  yield  the 
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Fig.  2.8.  Changes  in  intramolecular  bond  distances  with  pressure  at  300  K.  The  col¬ 
ors  correspond  to  different  compression  factors  (effective  pressures)  of  the  simulation 
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Fig.  2.9.  (a)  snapshots  from  ab  initio  MD  simulation  of  highly  compressed  liquid 
nitromethane  showing  a  proton  transfer  event,  (b)  Schematic  of  the  proton  transfer 
reaction  observed 


final  products  CH2N02~  and  CH3NO(OH)+.  The  H  transfer  event  occurs  on 
a  time  scale  of  about  15  fs. 

Since  the  C  N  bond  is  the  weakest  covalent  bond  in  nitromethane,  one 
might  think  that  the  C  N  bond  should  break  as  is  commonly  observed  in  the 
collision  simulations  of  the  previous  sections.  Indeed,  Fig.  2.8  reveals  that  as 
the  compression  of  the  liquid  is  increased,  there  is  a  significant  compression 
of  the  C— N  bond,  as  compared  to  the  CH  and  CO  bonds.  So  one  question 
to  ask  is  why  then  does  the  C  H  bond  first  break?  To  answer  this  question, 
we  used  molecular  orbital  theory  and  examined  the  electronic  structure  of  ni¬ 
tromethane.  As  previously  mentioned,  it  is  notable  that  there  is  a  significant 
compression  of  the  CN  bond  as  the  pressure  is  increased.  Examination  of  the 
frontier  molecular  orbitals  of  the  nitromethane  molecule  as  the  C  N  bond  is 
shortened  reveals  that  the  methyl  group  proton  of  nitromethane  becomes  a 
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Fig.  2.10.  Frontier  molecular  orbitals  of  nitromethane  as  a  function  of  the  CN  bond 
length.  The  top  portion  of  the  figure  show  the  calculated  orbital  energies  of  the 
highest  occupied  molecular  orbital  (HOMO)  and  the  lowest  unoccupied  molecular 
orbital  (LUMO).  The  bottom  portion  of  the  figure  shows  isosurface  plots  of  the 
LUMO  orbital  of  nitromethane.  As  the  CN  bond  is  compressed  the  LUMO  orbital 
energy  is  stabilized,  and  the  orbital  lobe  on  one  of  the  hydrogen  atoms  becomes 
larger,  as  highlighted  by  the  red  arrow 

better  Lewis  acid.  At  the  same  time,  the  carbonyl  oxygen  becomes  a  better 
Lewis  base.  This  is  illustrated  in  Fig.  2.10  which  shows  the  highest  occupied 
molecular  orbital  (HOMO)  and  lowest  unoccupied  molecular  orbital  (LUMO) 
energies  as  a  function  of  the  C  N  bond  distance.  The  HOMO  orbital  of  ni- 
tromethane  (not  shown)  can  be  described  as  a  pi-bonding  orbital  of  the  aro¬ 
matic  ONO  fragment  of  nitromethane.  This  orbital  is  also  polarized  toward  tin* 
oxygen,  and  as  the  C  N  bond  is  compressed,  it  increases  in  energy,  thereby 
making  the  oxygen  atoms  better  Lewis  base  sites.  The  LUMO  of  nitromethane. 
shown  on  the  bottom  of  Fig.  2.10,  can  be  described  as  an  anti-bonding  orbital 
of  the  aromatic  ONO  fragment.  However,  there  is  a  significant  component 
of  this  molecular  orbital  localized  on  the  hydrogen  atoms  as  highlighted  by 
the  red  arrow  in  Fig.  2.10.  Moreover,  this  contribution  becomes  larger  and 
the  LUMO  energy  becomes  stabilized  as  the  CN  bond  is  compressed.  Both  of 
these  show  that  the  nitromethane  protons  become  more  Lewis  acidic  as  the 
CN  bond  is  compressed. 

From  first  principles  molecular  dynamics  simulations,  Manaa  et  al.  [59] 
have  also  observed  proton  transfer  as  the  first  stage  of  the  nitromethane 
decomposition  reaction  at  high  temperature  (3,000  K)  and  high  density 
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(p  =  1 .97 g  cm-3)  conditions.  The  same  observation  was  made  under  shock 
compression  conditions  with  the  shock  speed  of  7  km  s_1  [60].  Some  early, 
static  high  pressure  experiments  [61]  show  that  the  thermal  detonation 
time  for  deuterated  nitromethane  is  ten  times  longer  than  that  of  ordinary 
nitromethane.  While  more  recent  experiments  comparing  the  detonation 
of  nitromethane  and  deuterated  nitromethane  show  less  dramatic  differ¬ 
ences  [56,62],  the  measured  particle  velocities  in  the  reaction  zone  of  deto¬ 
nation  reveal  that  changes  upon  deuteration  are  most  distinct  in  the  early 
part  of  the  reaction  zone.  Consistent  with  these  experimental  results,  our 
calculations  do  suggest  that  proton  transfer  might  be  important  in  the  initial 
stages  of  the  reaction.  Moreover,  the  quantum  chemical  analysis  of  the  proton 
transfer  step  as  a  function  of  the  CN  bond  distance  provides  fundamental 
insight  into  why  this  process  occurs  more  readily  as  the  pressure  is  increased. 


2.4  High  Pressure  Non-molecular  Solid  Phases 
of  Polynitrogen 

A  covalent,  single-bonded  nitrogen  phase,  analogous  to  carbon  in  diamond, 
has  long  been  thought  to  exist  at  high  pressures.  If  such  a  material  were 
metastable  at  ambient  pressures,  it  would  make  an  extremely  powerful  high- 
energy-density  material  (HEDM).  Due  to  the  uniquely  large  difference  in 
energy  between  the  N  N  single  and  triple  bonds,  the  energy  density  of  single 
bonded  polymeric  nitrogen  has  been  calculated  [63]  to  be  at  least  0.4  eV  cm-3, 
which  is  about  three  times  that  of  the  most  powerful  energetic  materials 
known  today. 

The  structure  of  single-bonded  polymeric  nitrogen  was  at  first  suggested 
to  be  that  of  the  known  single-bonded  forms  of  nitrogen’s  group  15  con¬ 
geners,  namely,  the  black  phosphorus  (BP)  structure  and  the  A7  structure 
of  arsenic.  In  1992,  based  on  possessing  gauche  interactions  among  the  lone 
pairs  of  adjacent  nitrogen  atoms,  Mailhiot,  Yang  and  McMahan  [64]  predicted 
a  unique  single-bonded  structure  that  has  no  analogue  in  other  natural  st  ruc¬ 
tures  that  they  called  the  cubic  gauche  ( CG)  structure.  Based  on  first  prin¬ 
ciples  calculations,  within  a  large  range  of  atomic  volumes,  the  CG  structure 
was  predicted  to  be  more  stable  than  the  ZIP  and  A  7  phases.  Further  computa¬ 
tional  studies  [65]  showed  that  CG  might  be  metastable  at  ambient  pressures, 
thereby  furthering  the  prospect  of  using  polymeric  nitrogen  ns  a  practical 
HEDM. 

Experimentally,  extensive  efforts  were  made  to  produce  polymeric  nitrogen 
and  evidence  for  non-molecular  amorphous  phases  of  nitrogen  at  HMKH’a 
and  1,000 K  were  reported  [66  68].  In  one  such  study  |6(i|,  the  amorphous 
material  was  even  recovered  at  ambient  pressures  and  low  temperature.  It  is 
likely  that  the  amorphous  materials  observed  were  mixtures  of  small  clusters 
of  non-molecular  phases.  In  2004,  Erements  and  coworkers  reported  [69]  the 
synthesis  of  a  crystalline  form  of  single- bonded  polymeric  nitrogen  at  pressures 
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greater  than  llOGPa  and  2,000  K  temperatures.  Under  these  conditions,  the 
molecular  nitrogen  is  transformed  into  the  cubic  gauche  crystal  phase.  The 
apparent  high  energy  barrier  between  molecular  and  non-molecular  phases 
was  overcome  by  the  high  temperature  synthesis.  A  strong  pressure-dependent 
hysteresis  was  observed  during  the  synthesis  [66,69].  The  high  sensitivity  of 
the  final  phase  on  its  pre-history  was  also  noted  in  first-principles  simulations 
[64,  70—72],  and  suggests  the  existence  of  many  competing  metastable  high 
pressure  phases  separated  by  large  energy  barriers. 

Recently,  a  number  of  new  phases  of  polymeric  nitrogen  have  been  found 
by  first  principles  molecular  dynamics  calculations  at  high  pressures  and  tem¬ 
peratures  [70-72].  Alemany  and  Martins  [71]  found  a  metastable  metallic 
phase  composed  of  zig-zag  chains  of  nitrogen  atoms  packed  in  a  body-centered 
orthorhombic  structure.  Through  similar  MD  simulations  using  temperatures 
(up  to  10,000  K),  Martin  and  coworkers  [70]  also  discovered  a  metastable  non- 
molecular  phase  composed  of  zig-zag  chains  of  alternating  N  N  single  and 
double  bonds.  Although  not  single- bonded,  this  new  zig-zag  chain  phase  is 
lower  in  energy  than  CG  at  low  pressures  and  competitive  with  CG  at  higher 
pressures. 

2.4.1  Polynitrogen  Phases  from  Simple  Cubic  Motifs 

Although  successful,  the  search  for  new  phases  by  first  principles  MD  sim¬ 
ulations  is  ad  hoc  in  nature  and  can  be  very  computationally  demanding. 
A  systematic  method  to  finding  new  polymeric  nitrogen  phases  is  obviously 
desirable.  The  fact  that  all  single-bonded  metastable  phases  of  group  15 
elements  predicted  or  found  in  nature  can  be  considered  Peierls-like  distor¬ 
tions  [73]  of  the'  simple  cubic  (SC)  structure  hints  at  such  a  method.  We  have 
recently  developed  a  systematic  method  of  looking  for  new  structures  of  solid, 
single-bonded  polymeric  nitrogen  [74].  The  method  is  a  synthesis  of  modern 
first  principle's  calculations  with  the  geometrical  model  of  crystal  structure 
developed  by  We  lls  [75  77]  and  Burdett  anel  McLarnan  [78].  The  procedure 
effectively  exploit*  the'  connection  of  the  other  structures  to  the  SC  reference 
structure,  allowing  focus  on  the  most  physically  meaningful  structure  candi¬ 
dates.  The'  procedure  is  able  to  recover  all  the  threefold  connected  nitrogen 
allotropes  found  before,  (e*.g.,  A7,  BP,  CG)  and  eight  new  metastable  single- 
bonded  allotropes. 

The  approach  presented  here  is  based  em  Peierls-like  distortions  erf  the 
SC  reference  structure  in  order  to  systematically  generate  new  single-bonded 
polynitrogen  structures  consistent  with  a  threefold  connectivity  pattern 
of  nitrogen.  In  the  undistorted  SC  structure,  each  nitrogen  atom  has  six 
equidistant  nearest  neighbors.  The  distortion  must  be  such  that  each  nitrogen 
atom  forms  only  thns*  covalent  bonds,  thereby  breaking  the  remaining  three 
of  the  ‘bonds’  of  the  original  5C1  structure.  To  obtain  a  complete  set  of  possible 
SC  distortions  consistent  with  threefold  connectivity  we  use  the  combinato¬ 
rial  analysis  procedure  of  Burdett  and  McLarnan  [78],  which  is  based  on  a 
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#1  -  BP  #3  -  CG  #1 1 

Fig.  2.11.  Connectivity  diagrams  (or  structure  types)  of  a  simple  cubic  reference 
structure  consistent  with  a  threefold  coordination  of  nitrogen.  Shown  are  the  con¬ 
nectivity  diagrams  of  structure  #1,  #3  and  #11  of  a  possible  36.  The  structure 
numbering  that  is  adopted  is  that  of  [78].  A  single  eight  atom  SC  unit  cell  is  high¬ 
lighted  with  red  atoms  and  bonds 


generalization  of  Polya’s  enumeration  theorem  [79,80].  In  particular,  Burdett 
and  McLarnan  show  that  out  of  4,096  possible  combinations  there  are  only  36 
unique  ways  of  labeling  the  edges  of  an  eight  atom  5Cunit  cell  as  “bonds”  and 
“non-bonds”  consistent  with  the  threefold  connectivity.  Those  bond-breaking 
patterns  can  be  described  by  connectivity  diagrams  or  structure  types,  as 
depicted  in  Fig.  2.11.  The  complete  specification  of  all  the  36  structure  types 
is  provided  in  [78]  and  we  adhere  to  the  same  numbering.  Note  that  all  the 
previously  reported  structures  of  single-bonded  nitrogen,  such  as  CG,  BP, 
A7,  and  LB  [72]  are  covered  in  the  set  of  36  structure  types.  Figure  2.11 
shows  the  connectivity  diagrams  for  BP  (structure  #1),  CG  (#3)  and  a  new 
structure  #11  that  will  be  discussed  in  further  detail  later. 

Once  a  structure  type  is  chosen,  the  next  step  is  to  use  the  connectivity 
pattern  to  distort  the  SC  structure  to  generate  an  initial  geometry  that  can  be 
optimized  by  first  principles  calculations.  The  most  straightforward  distortion 
is  to  displace  each  atom  of  the  unit  cell  by  a  distance,  d ,  along  the  body 
diagonal  toward  the  three  connected  neighbors  as  shown  in  Fig.  2.12.  Once 
all  atoms  are  distorted  in  this  way,  a  rough  initial  structure  is  generated  as 
depicted  on  the  r.h.s.  of  Fig.  2.12.  This  can  be  described  as  a  Peierls-like 
distortion  of  the  SC  structure.  This  initial  structure  acts  as  a  starting  point 
for  a  series  of  first  principles  calculations.  In  our  calculations,  we  utilized  the 
Kohn  Sham  density  functional  theory  calculations  with  the  Perdew  Burke- 
Ernzerhof  (PBE)  [27]  exchange-correlation  functional.  The  Vienna  ab  initio 
simulation  package  (VASP)  [81]  was  used  with  the  projector  augmented  wave 
(PAW)  method  of  Blochl  [82]  to  treat  the  core  states.  A  plane-wave  cutoff  of 
39  Ry  was  used  and  Brillouin-zone  integration  was  performed  using  a  8  x  8  x  8 
Monkhorst  Pack  gamma  point  centered  grid. 

Although  the  generation  of  the  initial  structures  may  seem  straightfor¬ 
ward,  the  structures  may  not  be  stable  at  all  pressures,  if  at  all.  Thus,  a  careful 
structure  optimization  procedure  is  needed  in  order  to  find  the  appropriate 
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Fig.  2.12.  Schematic  illustration  of  the  procedure.  Left:  a  structure  type  of  an  eight 
atom  cube  is  first  chosen  (blue  atoms  with  red  bonds).  A  Peierls-like  distortion  that 
is  consistent  with  the  connectivity  is  applied.  Right:  the  structure  after  distortion 
of  all  atoms  in  an  extended  representation 


ranges  of  pressure/ volume  conditions,  as  well  as  to  prevent  premature  decom¬ 
position  of  structures  before  they  get  sufficiently  close  to  the  final  optimized 
geometries. 

Here  we  adopted  the  following  computational  procedure.  Beginning  with 
the  initial  structure  generated  by  the  aforementioned  Peierls-like  distortion, 
the  lattice  constant  a  is  scanned  from  4.1  A  downwards,  with  a  step  of  0.05  A. 
For  each  value  of  a,  the  following  optimization  steps  are  performed. 

•  First  the  distortion  parameter  d  is  adjusted  to  minimize  the  ab  initio 
calculated  total  energy. 

•  Then  the  unit  cell  shape  is  allowed  to  relax,  fixing  the  fractional  coordi¬ 
nates  of  the  atoms  and  the  unit  cell  volume. 

•  Finally,  full  relaxation  of  both  the  unit  cell  shape  and  the  atomic  positions 
at  fixed  unit  cell  volume  a3  are  performed. 

At  each  such  relaxation  stage  the  structure  must  still  preserve  the  intended 
connectivity  pattern.  The  scan  of  the  lattice  constant  continues  until  either 
the  first  time  a  fully  optimized  structure  with  a  required  connectivity  pat¬ 
tern  is  encountered  or  a  =  3.0  A  is  reached.  In  cases  where  a  fully  optimized 
structure  that  preserves  the  intended  structure  type  is  obtained,  which  hap¬ 
pened  for  26  out  of  36  types,  mechanical  stability  of  the  structure  is  verified 
by  calculating  the  Phonon  spectrum,  and  also  the  enthalpy- pressure  curve. 
The  mechanical  stability  is  inferred  from  phonon  densities  of  states  for  each 
of  the  26  structures  at  various  pressures.  The  stabilities  are  deduced  from  a 
force  constant  matrix  calculated  in  a  2  x  2  x  2  (64  atom)  supercell  by  means 
of  consecutive  displacements  of  each  atom  of  the  original  unit  cell  by  0.02  A 
along  each  axis.  In  some  cases  this  supercell  size  may  be  inadequate,  giving 
rise  to  false  negative  frequencies.  In  this  respect,  only  the  affirmative  outcomes 
of  the  metastability  tests  should  be  treated  as  conclusive. 
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Fig.  2.13.  Results  of  the  phonon-spectra  based  mechanical  stability  tests  for  tin*  :t(i 
structure  types  and  different  pressures.  Mechanically  stable  and  unstable  structures 
at  different  pressures  are  shown  with  circles  and  crosses,  respectively  [74] 


Analysis  of  the  mechanical  stability,  shows  that  at  least  8  of  the  2G  new 
structures  are  mechanically  stable  at  pressures  less  than  100  GPa  (#6,  8,  9,  1 1 , 
1G,  20,  21  and  30).  Structures  #11  and  #30  are  even  found  to  be  mechanically 
stable  at  zero-pressure.  The  tests  of  mechanical  stability  at  various  pressures 
are  summarized  in  Fig.  2.13.  Figure  2.14  reveals  that  the  cubic  gauche  phase 
is  still  the  lowest  enthalpy  structure  over  the  whole  pressure  range  studied. 
Moreover,  it  also  shows  pronounced  mechanical  stability  in  our  phonon  cal¬ 
culations,  which  is  consistent  with  the  same  type  of  analysis  in  [65].  Leaving 
the  cubic  gauche  case  aside,  the  data  in  Fig.  2.13  suggest  that  new  structures 
generated  by  Peierls  distortions,  such  as  structure  types  #8  and  #11,  have 
(  luigies  that  are  comparable  or  lower  than  those  of  the  previously  studied 
stint  turcs  #1  ( BP ),  #7  {AT),  and  #15  {LB).  The  3-dimensional  structures 
of  structures  #8  and  #11  are  depicted  in  Fig.  2.15. 

To  get  an  insight  into  the  relative  energy  rankings  of  different  struct  ure 
type's  we  analyze  a  correlation  between  the  first-principles  total  energies  ami 
t  he  number  of  cis-,  gauche -,  trails-  dihedral  angles,  and  four-member  rings 
of  those  structure  types.  Based  on  such  analysis,  the  lowest  energy  structures 
of  the  low  pressure  region,  structures  #3  ( CG )  and  #11  (new)  are  recognized 
as  structures  that  minimize  the  number  of  /,ran.v-dihedrals  and  fdiii  membei 
lings  while  maximizing  the  //aur7m-dihedrals.  It  is  interesting  to  note  that 
the  stronger  weight  of  trans-  than  the  cis- dihedrals  is  in  contrast  with  the 
simplified  “lone-pair”  orbital  picture  [64,78]  gained  from  ah  init  io  ralculat  ions 
on  small  nitrogen  molecules. 

After  mechanical  stability  screening,  newly  found  struct iires  //  I  I  and  //8 
are  the  second  lowest  enthalpy  structures  after  cubic  gam  In  at  low  ( •  l()(  ;pa) 
and  nioderatc/high  (~>40GPa)  pressures,  respectively.  Selected  views  of 
structures  #11  and  #8  are  shown  in  Fig.  2.16.  Structure  //II  is  distinguished 
by  the  high  R,T2  symmetry  and  its  structure  can  be  described  by  sets  of 
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Structure  type  # 


Fig.  2.14.  Calculated  (a)  ent  halpy,  (b)  total  energy,  and  (c)  volume  per  atom  ver¬ 
sus  structure  type  #  at  pressures  of  lOGPa  (circles)  and  100  GPa  (crosses).  The 
enthalpies  and  energies  are  with  respect  to  those  of  the  CG  phase  at  the  corre¬ 
sponding  pressures.  Note  that  structures  #1,  3,  7,  and  15  are  the  previously  known 
phases  Bl\  GG,  A7,  and  LB ,  respectively  [74] 


(a)  (b) 

Fig.  2.15.  Kx tended  structures  of  two  now  low  energy  and  metastable  phases*, 
(a)  structure  #11  and  (b)  structure  #8 

cis  trans  chains  directed  along  the  three  coordinate  axes.  In  fact  then'  are 
two  mutually  orthogonal  sets  of  cis  trans  chains  along  each  axis.  Structure 
#8  lias  P_i  symmetry  and  a  very  strong  structural  resemblance  to  structure 
#11,  being  similarly  shaped  by  sets  of  cis  trans  chains  but  along  two  out  of 


2  First  Principles  Molecular  Simulation  of  Energetic  Materials 


95 


Fig.  2.16.  A  perspective  view  of  the  lattice  structure  of  the  Layered-Boat  (Hi) 
structure 


the  three  axes,  hence  its  lower  symmetry.  Structure  #8  contains  large  voids  in 
between  those  layers  which  may  be  responsible  for  the  less  favorable  enthalpy 
at  lower  pressure.  Structure  #1 1  is  an  insulator  with  an  indirect  band  gap  de¬ 
creasing  from  2.8 eV  at  near  zero  pressure  to  2.2 eV  at  43  GPa.  Structure  #8  is 
also  an  insulator  with  an  indirect  gap  at  low  pressures  (e.g.,  l.OeV  at  26GPa) 
replaced  by  a  higher  energy  direct  gap  at  higher  pressures  (1.5 eV  at  96GPa). 


2.4.2  Polynitrogen  Phases  from  Chain  Motifs 

Using  a  series  of  high  temperature  and  high  pressure  ab  initio  molecular  dy¬ 
namics  simulations  we  have  previously  found  [72]  a  new  single  bonded  non- 
molecular  phase  that  we  have  named  the  ‘layered  boat.’  or  LB  phase  as  shown 
in  Fig.  2.16.  The  new  phase  consists  of  layers  of  fused  six-membered  rings  in  a 
boat  conformation,  and  thus  the  layered  boat  name.  LB  has  a  P2i/m  symme¬ 
try  and  consists  of  two-dimensional  layers  of  fused  N^-rings  where  the  rings 
are  all  in  the  boat  conformation.  A  calculated  enthalpy  versus  pressure  phase 
diagram  reveals  that  LB  lies  between  BP  and  A  7  in  enthalpy  at  pressures  less 
than  200  GPa.  At  pressures  higher  than  this,  the  enthalpy  of  LB  edges  above 
that  of  A  7. 

The  structure  of  LB ,  initially  obtained  at  high  pressures,  remains  quali¬ 
tatively  unchanged  when  re-optimized  at  different  pressures  in  the  range  of 
0  300 GPa,  indicating  that,  along  with  other  polymeric  phases  of  nitrogen, 
LB  is  also  metastable  for  a  wide  range  of  pressures.  An  enthalpy  versus  pres¬ 
sure  phase  diagram  for  LB  and  various  other  phases  is  shown  on  f  ig.  2.17 
for  pressures  less  than  200 GPa.  At  low  pressures,  from  OGI’A  to  50 GPa, 
zzCH  is  predicted  to  be  the  lowest  in  enthalpy,  in  agreement  with  I  lie  work  of 
Mattson  et  al.  [70].  At  pressures  greater  than  50 GPa  and  less  than  200 GPa, 
CG  becomes  the  lowest  in  enthalpy,  and  BP  becomes  lowest  at  even  higher 
pressures.  The  new  phase,  LB ,  is  higher  in  enthalpy  than  ( 'G  and  zz CH  for 
all  observed  pressures  and  is  between  A  7  and  BP  for  pressures  up  to  210  GPa. 
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Fig.  2.17.  Phase  diagram  of  various  polymeric  nitrogen  phases.  The  enthalpy  per 
atom  is  relative  to  the  enthalpy  per  atom  of  CG  at  100  GPa.  For  clarity,  the  phase 
diagram  is  split  into  two  pressure  regions:  (a)  0-80 GPa  and  (b)  80-200 GPa  [72] 


At  pressures  above  210 GPa,  LB  becomes  slightly  higher  in  enthalpy  than  AT. 
SC  is  the  least  favorable  phase  among  all  the  phases  discussed  here.  CH  has 
the  fastest  growth  of  enthalpy  with  respect  to  the  increasing  pressure,  and 
quickly  becomes  less  favorable  than  all  except  SC  as  the  pressure  is  increased. 

The  LB  phase  can  be  constructed  using  the  SC  starting  structures  as 
described  in  the  previous  section.  Viewed  from  an  alternative  perspective,  the 
structure  of  the  LB  phase  is  related  to  the  BP  and  A  7  phases  in  that  they  are 
also  composed  of  layers  of  fused  six-member  rings,  however,  in  the  latter  two 
the  rings  are  in  the  chair  conformation.  The  LB,  BP  and  AT  phases  can  be 
thought  of  as  different  arrangements  of  the  zig-zag  chains  in  the  zzCH  phase 
approaching  along  given  directions  and  forming  new  inter-chain  bonds.  This  is 
shown  in  Fig.  2.18  where  the  zzCH  phase  is  depicted  with  possible  inter-chain 
bonds  that  would  result  in  the  formation  of  the  LB,  BP  and  A  T  phases.  The 
two-dimensional  layers  emerge  with  a  periodic  six-ring  pattern  as  a  result  of 
the  chains  in  zzCH  coming  close  together  and  forming  bonds  which  change 
the  coordination  of  the  nitrogen  atoms  from  two  to  three. 

In  this  representation,  the  ‘chains’  of  A  T  art*  aligned  parallel  to  each  other, 
and  all  are  simultaneously  tilted  to  one  side.  In  Fig.  2.18  the  chains  arc* 
depicted  tilted  to  the  left.  An  enantiomeric  ‘right  Lilted’  form  of  A  7  could 
be  constructed  from  zzCH  by  choosing  a  different  direction  along  which  the 
zig-zag  chains  come  together.  In  the  same  spirit,  LB  could  be  characterized 
as  an  arrangement  of  zig-zag  chains  tilted  to  the  left  and  right  in  an  alternat¬ 
ing  fashion.  The  relationship  between  the  layered  phases  and  the  zzCH  phase 
outlined  here  is  not  meant  to  show  how  the  layered  phases  are  most  favor¬ 
ably  formed  from  the  zzCH  phase.  Nonetheless,  it  provides  possible  pathways 
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Fig.  2.18.  The  zzCH  phase  with  the  possible  inter-chain  bonding  drawn:  BP  (dash- 
dot  lines),  A  7  (dash  lines)  and  LB  (dot  lines) 


gauche  trans  cis 


Fig.  2.19.  Newman  projections  of  the  RjN  NRj  fragment  in  various  conforma¬ 
tions.  ‘lp’  represents  the  lone-pair 


for  the  formation  of  the  layered  phases  from  the  zzCH  phase  that  do  not  re¬ 
quire  significant  structural  rearrangement.  We  note  that  both  BP  and  LB  can 
also  be  considered  in  terms  of  cis  trans  chains  of  the  CH  phase  coming  close 
together  and  forming  inter-chain  bonds. 

The  three'  layered  phases,  A7,  BP,  and  LB  all  consist  of  fused  six-member 
rings  composed  of  nitrogen  atoms  that  are  spA  hybridized.  The  local  structure 
of  the  three  phases  is  distinct,  for  example,  the  rings  in  .A 7 and  BP  are  in  the 
chair  conformation  whereas  they  are  in  the  higher  energy  boat  conformation 
in  LB.  Here  we  examine'  the  relationship  between  the  local  structure  and 
the  relative  internal  energies  of  the  phases  by  considering  t  he*  dihedral  angles 
between  N  N  single  bonds. 

Quantum  chemical  calculations  of  small  molecules  of  t  he  form  IQN  NIQ 
(R  =  H,  CH:t,  etc)  show  that  they  have  two  minimum  energy  con  formers 
gauche  and  trans  as  depicted  in  Fig.  2.19.  The  gauche  struct  me  of  IQN  NIQ 
is  not  ideal  in  that  the  lp  N  N  lp  dihedral  angle,  9,  is  not  00°  but  90°. 
The  90°  dihedral  angle  of  the  gauche  conformer  can  be  accounted  for  by  the 
fact  that  this  0  angle  minimizes  the  two-orbit al/ four-electron  destabilizing 
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interaction  between  the  two  lone  pair  orbitals  of  the  adjacent  N  atoms. 
(We  note  that  gauche  definition  sometimes  used  is  somewhat  different  than 
the  convention  used  here.  For  example,  the  lp  N  N  lp  dihedral  angle,  in 
the  cubic  gauche  phase  is  107°  with  a  lp  N  N  R  dihedral  angle  being  close 
to  zero.)  For  a  later  discussion  we  also  introduce  the  cis  conformation  where 
the  lp — N  N  lp  =  0°.  This  is  actually  the  maximum  energy  conformation 
as  the  lp  N  N  lp  is  rotated  through  360°. 

The  energy  of  the  N — N  moiety  depends  on  the  mutual  orientation  of 
two  adjacent  lone-pairs  as  defined  by  the  lp  N  N  lp  dihedral  angle.  The 
structure  of  CG  was  postulated  [64]  as  a  solid  phase  of  nitrogen  that  has  all  of 
its  single  bonds  in  a  gauche  dihedral  conformation.  Indeed,  our  calculations 
as  well  as  many  others  [64,70,71,83]  demonstrate  that  above  50  GP,  CG  has 
the  lowest  enthalpy  among  all  the  known  nitrogen  phases. 

The  internal  energy  of  the  A7,  BP ,  and  LB  phases  can  be  related  to  the 
ratio  of  gauche. ,  tiuns  and  cis  N  N  bonds  there  are  in  each  pluus<\  There  arc* 
six  such  dihedral  angles  in  each  six-member  ring  of  these  phases.  BP  has  a 
ratio  of  4-gauche/2-trans/0-cis,  while  .4  7  has  a  ratio  of  W gauche/ 6-  trans/iY  cis 
and  LB  has  a  ratio  0-gauche/  4-trans/2-cis.  With  the  gauche  conformation 
being  the  most  stable,  followed  by  the  trans-c onformation  and  the  cis  confor¬ 
mation  being  the  least  stable,  it  is  clear  that  this  simple  structural  analysis 
predicts  the  following  ordering  of  the  three  layered  six-member  ring  phases: 
BP  <  Al  ^  L JB .  First- principles  computations  of  the  three  phases  aie  m 
agreement  with  this.  The  calculated  internal  energy  per  atom  of  the  phases 
is  —268.6710 eV  <  -268.6103eV  <  -268.5912eV  at  20GPa.  Of  course,  the 
overall  enthalpy  of  the  phases  contains  a  volume  factor,  which  has  substantial 
contribution  at  elevated  pressures. 

2.4.3  Polynitrogen  Phases  from  Helical  Motifs 

Helical  chains  can  also  be  considered  as  a  natural  structural  theme  in  poly¬ 
meric  nitrogen  phases.  Recognition  of  helical  structural  motifs  in  the  exper¬ 
imentally  observed  CG  crystal  lattice  has  lead  to  the  discovery  of  a  new 
single-bonded  non- layered  nitrogen  structure  named  chairvd  web  (C'W)  [84]. 
First-principles  density  functional  theory  calculations  reveal  that  GW,  which 
was  originally  identified  at  high  pressures  is  metastable  at  ambient  conditions 
as  well.  The  metastability  is  demonstrated  by  both  high-quality  phonon  dis¬ 
persion  calculations  and  finite  temperature  first-principles  molecular  dynamics 
simulations.  In  addition  the  new  GIF  phase  is  thermodynamically  more  stable 
than  the  CG  phase  in  the*  ambient  pressure  regime. 

The  results  were  obtained  by  first  principle's  density  functional  theory  cal¬ 
culations  with  the  Perdew  Burke  Ernzerhof  (PBE)  [27]  exchange-correlation 
functional.  The  SIESTA  [85]  package  was  used  for  most  of  the*  exploratory 
simulations,  preliminary  optimization  of  structure's,  and  molecular  dynam¬ 
ics  simulations,  while  VASP  [81]  was  used  to  calculate*  the  final  phase*  dia¬ 
gram.  phonon  dispersion  and  band  structure.  For  the  SIESTA  calculations, 
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the  Troullier-Martins  [86]  norm-conserving  pseudopotential  with  a  nitrogen 
reference  configuration  of  [He]2s22p*  and  a  cut-off  radius  of  0.98  A  was  uti¬ 
lized.  A  custom  numerical  4double-zeta  with  polarization’  SIESTA-type  basis 
set  [87]  was  developed  for  the  calculations.  A  10  A  cut-off  was  employed 
for  the  fc-point  sampling.  The  calculations  performed  with  VASP  were  based 
on  the  projector  augmented  wave  (PAW)  method  of  Blochl  [82]  where  a 
39  Ry  plane  wave  cutoff  and  Brillouin-zone  integration  with  12  x  12  x  12 
Monkhorst-Pack  grid  were  used.  In  both  programs  a  variable-cell-shape  con¬ 
jugate  gradient  method  under  constant  pressure  was  used  for  the  minimization 
of  the  geometries,  and  molecular  dynamics  was  performed  within  the  Nose- 
Parinello-Raman  scheme. 

Ab-initio  calculations  on  small  Molecules  of  the  type  R2  N  N  R2 
demonstrated  that  the  most  favorable  dihedral  angle  is  the  so-called  gauche 
dihedral,  i.e.,  ~120°,  due  to  the  effect  of  hyperconjugation  and  mutual  repul¬ 
sion  of  nitrogen  lone  electron  pairs  [64] .  A  chain  of  atoms  constructed  with 
all  gauche  dihedral  angles  results  in  the  formation  of  right-  or  left-handed  he¬ 
lices,  depending  on  whether  a  dihedral  of  -(-120°  or  -120°  is  applied.  Indeed, 
within  the  CG  structure,  helices  that  have  four  and  eight  atoms  per  turn  can 
be  identified  along  the  [100]  direction  and  threefold  helices  can  be  located 
along  the  [110]  direction.  Figure  2.20d  highlights  the  eightfold  helices  in  CG. 

Pure  helices  are  also  observed  in  some  high  pressure  phases  of  group 
16  elements,  such  as  sulfur,  selenium  and  tellurium  [88-90]  as  well  as  in 
scandium  [91].  Unlike  the  group  16  elements,  nitrogen  cannot  form  stable 
helices  on  its  own.  In  order  to  acquire  stability  a  single  helical  chain  of  nitro¬ 
gen  atoms  has  to  have  its  valence  saturated  (form  three  single  bonds).  If  the 
valence  is  not  saturated,  then  partial  double  bonds  will  form  and  the  chain 
will  flatten  out  to  either  a  cis  trans  or  zig-zag  chain.  In  the  latter  case  the 


Fig.  2.20.  (a)  Flipping  or  ‘inverting’  two  (highlighted  in  red)  out  of  four  N2  mole¬ 
cules  in  the  primitive  unit  cell  of  the  molecular  £-phase  (l.h.s.)  results  in  an  eightfold 
helical  chain  structure  (r.h.s.)  (b)  £-phase,  (c)  ‘inverted’  £-phase  and  (d)  CG.  For  B, 
C  and  D  a  2  x  2  x  2  supercell  is  shown  with  one  helical  chain  highlighted  in  red  [84] 
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energy  lowers  by  0.5  eV  per  atom.  The  reverse  tendency  also  holds,  namely 
cis  trans  or  zig-zag  chains  acquire  helicity  once  capped  due  to  the  above  men¬ 
tioned  'gauche-effect’  and  the  energy  is  lowered  by  0.6 eV  per  capped  atom.  It 
is  interesting  to  note  that  in  both  cases  the  lowering  of  energy  is  achieved  by 
changes  in  the  local  hybridization  (from  sp 3  to  sp2  and  vice  versa),  symmetry 
breaking  and  splitting  of  degeneracy  at  the  Fermi  level.  These  are  characteris¬ 
tic  of  a  Peierls  deformation.  If  a  pure  nitrogen  structure  is  sought,  the  helices 
have  to  be  brought  close  enough  to  each  in  order  to  create  inter-helical  bonds. 
The  geometry  of  CG  could  be  thought  of  as  such  a  structure  where  interlinked 
nitrogen  helices  form  a  single-bonded  network  as  illustrated  in  Fig.  2.20d. 

The  base  helical  structure  of  CG  can  be  formed  from  what  has  been 
proposed  to  be  the  molecular  £-phase,  which  itself  is  considered  to  be  the 
immediate  molecular  precursor  [92]  of  CG  during  its  synthesis.  Helices  can 
be  formed  by  changing  the  alignment  of  two  out  of  the  four  N2  molecules 
within  the  primitive  unit  cell  of  the  (^-structure.  The  necessary  rotation  about 
the  two  molecule’s  center  of  mass  as  to  'invert’  their  orientation  is  shown  in 
Fig.  2.20a.  Using  this  ‘inverted’  (^-structure  (Fig.  2.20b)  as  a  starting  geome¬ 
try  for  a  first-principles  calculation,  eightfold  helical  chains  form  after  only  a 
few  optimization  steps  as  depicted  in  Fig.  2.20c.  Further  optimization  at  high 
pressure  leads  to  the  formation  of  the  single-bonded  CG  phase  (Fig.  2.20d). 
In  a  similar  manner  the  geometry  optimization  of  properly  arranged  four-  and 
threefold  helical  chains  also  result  in  the  formation  of  CG. 

It  is  natural  to  propose  that  certain  arrangements  of  the  helical  nitrogen 
chains  could  form  the  basis  of  other  polynitrogen  structures.  To  achieve  a 
structure  where  all  nitrogen  atoms  are  three-coordinate,  bonds  between  adja¬ 
cent  chains  must  be  made  because  the  nitrogen  atoms  are  only  two-coordinate 
in  the  pure  helical  chains.  Different  arrangements  of  3-,  4-,  etc.  fold  helices 
or  of  their  combinations  could  lead  to  alternative  structures.  For  example, 
a  slight  rearrangement  of  the  threefold  helices  within  CG  results  in  a  struc¬ 
ture  that  has  recently  been  proposed  as  a  new  polynitrogen  structure  [93]. 
Arrangement  and  connection  of  sixfold  helical  chains  as  shown  in  Fig.  2.21, 
followed  by  geometry  optimization,  produced  a  new  stable  structure  we  have 
called  chaired-web  ( CW ).  Connecting  the  helical  chains  together,  results  in  the 
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Fig.  2.21.  A  properly  arranged  set  of  sixfold  helices  (l.h.s.)  forming  the  single- 
bonded  CW  phase  (r.h.s)  [84] 
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(a)  (b) 

Fig.  2.22.  Lattice  structure  of  CW:  (a)  Perspective  view,  (b)  A  view  along  [111]. 
Left-  (red)  and  right-  (green)  oriented  sixfold  helices  with  six-rings  (white)  lying 
between  them  within  CW  phase  (3x3x3  supercell)  [84] 


formation  of  a  network  of  six-member  rings  that  are  in  the  so-called  ‘chair’ 
conformation,  as  the  name  given  to  the  structure  suggests.  Figure  2.21a  high¬ 
lights  the  six- member  rings  of  the  CW. 

The  structure  and  symmetry  of  the  new  phase  is  better  seen  in  Fig.  2.22. 
If  the  original  helices  turn  clockwise  (depicted  in  Fig.  2.22b  in  red)  then  new 
helices  with  counter-clockwise  turns  (depicted  in  Fig.  2.22b  in  green)  are  natu¬ 
rally  formed  in  between  the  original  helices  due  to  the  inter-helical  bonds.  The 
aforementioned  six-member  rings  are  formed  in  the  voids  between  the  helices 
as  shown  in  white  in  Fig.  2.22b.  These  voids  accommodate  the  nitrogen  lone 
pairs.  All  of  the  six-member  rings  are  in  the  so-called  chair  conformation.  The 
interplay  of  helical  chain  and  six-ring  structural  motifs  in  CW  is  reminiscent 
of  cis-trans  chain  and  six-member  ring  motifs  interplay  in  the  layered  phases 
DP A  7  and  LB  discussed  above  [72].  Additionally,  we  note  that  the  structure 
of  CW  can  be  obtained  by  applying  the  Peierls  distortion  method  for  generat¬ 
ing  purely  single-bonded  polymeric  nitrogen  [74].  However,  the  8-atom  simple 
cubic  reference  cell  first  demonstrated  [74]  would  not  be  large  enough  to  gen¬ 
erate  CW.  Rather  a  significantly  expanded  48  atom  (or  larger)  simple  cubic 
reference  structure  would  have  to  be  employed. 

The  primitive  unit  cell  of  CW  is  rhombohedral.  At  an  intermediate  pressure 
of  28GPa  the  unit  cell  vector  length  is  3.5  A  with  an  angle  measuring  99.172° 
and  the  corresponding  symmetry  group  is  R  3m.  The  six  equivalent  atomic 
positions  obtained  by  the  use  of  symmetry  in  fraction  coordinates  are  (0.8756, 
0.8756,  0.3206).  There  are  two  types  of  interatomic  distances  in  the  crystal 
lattice,  a  shorter  N  N  distance  within  a  given  six-member  ring  and  slightly 
longer  N  N  distance  between  two  adjacent  six-member  rings.  At  28GPa 
these  distances  are  1.36  A  and  1.45  A  respectively.  There  is  a  2:1:0  ratio  of 
gauche trans-  and  cis-  dihedral  angles  in  the  crystal  lattice.  Band  structure 
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Fig.  2.23.  Relative  enthalpy  versus  pressure  phase  diagram  of  CW ,  CG,  BP  and 
>17.  Enthalpy  is  relative  to  CG  at  0°  pressure.  The  inset  in  the  lower  right  corner 
enlarges  the  low-pressure  region  [84] 

calculations  show  that  CW  is  an  insulator  at  low  pressure  with  a  calculated 
band  gap  of  ~  5eV.  The  phase  diagram  shown  on  Fig.  2.23.  demonstrates  how 
close  CW  is  to  CG  in  terms  of  enthalpy.  At  zero  pressure  the  enthalpy  favors 
CW  over  CG  by  approximately  20meV.  The  temperature-dependant  free- 
energy  and  zero-point  energy  corrections  within  the  harmonic  approximation 
have  been  evaluated  at  zero  pressure.  The  zero-point  energy  correction  favors 
CG  over  CW  by  54  meV  thus  reversing  the  energy  order  of  the  two  phases  at 
zero  temperature.  With  the  raise  in  temperature  there  is  a  cross-over  of  the 
two  free-energy  dependencies  around  200  K.  At  ambient  temperature  the  new 
CW  phase  is  thermodynamically  more  stable  than  the  CG  phase  with  83meV 
in  free  energy. 

The  results  of  high-quality  phonon  dispersion  calculations  reveals  that  the 
CW  phase  is  metastable  at  very  low  pressures.  Although  the  lack  of  negative 
frequencies  in  the  phonon  dispersion  calculation  is  theoretically  rigorous  veri¬ 
fication  of  CW s  metastability  (at  least  at  low  temperatures),  three  additional 
tests  to  probe  the  energy  barrier  of  decomposition  have  also  been  performed: 

(i)  a  random  displacement  test  in  which  the  structure  is  optimized  after  the  po¬ 
sitions  of  atoms  are  randomly  displaced  by  5%  of  their  inter-atomic  distance; 

(ii)  ab-initio  molecular  dynamics  test  at  OGPa,  300  K  for  10  ps.  (iii)  a  hybrid 
test  in  which  ten  random  snapshots  of  the  initially  equilibrated  system  are 
taken,  the  positions  and  velocities  are  randomly  perturbed,  and  the  system  is 
further  a  subject  to  AIMD  at  OGPa,  300 K  for  an  additional  4ps.  A  rigorous 
study  of  CW s  metastability  lifetime  is  substantially  more  challenging,  but 
the  tests  indicate  that  the  lifetime  of  CW  may  be  of  practical  interest. 

The  structure  of  CW  is,  as  are  other  purely  single-bonded  polymeric 
nitrogen  phases,  governed  to  a  large  extent  by  the  repulsion  of  the  lone  pairs. 
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In  CW  the  inherent  stability  of  the  initial  gauche  helices  and  the  newly  formed 
helices  of  reverse  handedness,  is  offset  by  the  repulsion  of  the  lone  pairs  from 
neighboring  helices.  Importantly,  in  CW  the  lone  pairs  point  into  large  void 
spaces,  thereby  reducing  the  repulsion.  At  high  pressure  the  large  voids  are 
thermodynamically  unfavorable  because  of  the  associated  enthalpic  penalty. 
On  the  other  hand,  at  low  pressure,  the  balance  between  these  two  contradic¬ 
tory  influences  of  the  lone  pairs  is  such  that  CW  is  more  stable  than  CG  at 
ambient  pressures. 

In  looking  for  a  possible  precursor  of  CW  one  has  to  note  that  the  tra¬ 
ditional  understanding  of  N2  phase  diagram  was  recently  questioned  and 
preliminary  results  suggest  the  need  for  its  extension  and/or  revision  [94]. 
In  particular,  Gregoryanz  et  al.  [95]  found  a  molecular  phase  characterized 
by  a  configuration  of  three  nitrogen  molecules  nearly  forming  a  six-member 
ring.  One  could  speculate  that  CW  could  be  formed  from  such  a  phase  under 
appropriately  engineered  physical  conditions. 


2.5  Final  Remarks 

The  results  presented  in  this  chapter  indicate  that  significant  progress  can 
be  achieved  in  applying  various  first-principles  computational  methods  for 
atomic-scale  descriptions  of  energetic  systems  in  gas,  liquid  or  solid  phases. 
Using  well  developed  computational  chemistry  methods  it  has  been  possible 
to  determine  accurate  information  about  the  structural  properties  of  various 
systems,  the  relative  stabilities  of  different  configurations  in  different  phases, 
and  the  decomposition  reaction  mechanisms.  These  types  of  data  allowed  a 
direct  correlation  and  interpretation  of  the  corresponding  experimental  data. 
Moreover,  in  many  instances  the  type  of  data  obtained  theoretically  has  sub¬ 
stituted  for  the  lack  of  information  achievable  through  experimental  means. 
This  is  particularly  the  case  for  example,  for  description  of  complicated  reac¬ 
tion  mechanisms  where  accurate  identification  of  transition  states  for  various 
pathways  is  experimentally  very  challenging. 

Beside  structural  and  energetic  data,  important  steps  have  also  been 
achieved  computationally  in  description  of  the  dynamic  processes  in  condensed 
phases.  Among  the  computational  methods  used,  first  principles  calculation 
methods  are  the  first  choice  to  accurately  evaluate  the  structural,  energetic, 
electronic,  and  even  spectral  data. 

Alternatively,  the  progress  made  in  computational  hardware  mid  in  compu¬ 
tational  algorithms  should  allow  further  development  of  combined  quantum 
mechanical /molecular  mechanical  methods  for  description  of  energetic  and 
dynamical  problems,  for  systems  of  increased  size  and  complexity.  The  use 
of  ab  initio  molecular  dynamics  methods  for  description  of  the  properties  of 
ionic:  systems  is  another  area  which  is  expecting  to  grow  in  importance  and 
relevance. 
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